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Abstract. A useful concept for finding numerically the dominant correlations of 
a given ground state in an interacting quantum lattice system in an unbiased way 
is the correlation density matrix. For two disjoint, separated clusters, it is defined 
to be the density matrix of their union minus the direct product of their individual 
density matrices and contains all correlations between the two clusters. We show how 
to extract from the correlation density matrix a general overview of the correlations 
as well as detailed information on the operators carrying long-range correlations and 
the spatial dependence of their correlation functions. To determine the correlation 
density matrix, we calculate the ground state for a class of spinless extended Hubbard 
models using the density matrix renormalization group. This numerical method 
is based on matrix product states for which the correlation density matrix can be 
obtained straightforwardly. In an appendix, we give a detailed tutorial introduction 
to our variational matrix product state approach for ground state calculations for 1- 
dimensional quantum chain models. We show in detail how matrix product states 
overcome the problem of large Hilbert space dimensions in these models and describe 
all techniques which are needed for handling them in practice. 
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1. Introduction 



In an interacting quantum lattice model the ground state may have several kinds of 
correlations, such as long-range order, power-law, or exponentially decaying correlations. 
In the numerical treatment of such a model it is not clear a priori what kind 
of correlation will be dominant and what kind of operators corresponds to these 
correlations. Before calculating correlation functions, one typically chooses in advance 
which operators to consider, using prior knowledge and making initial assumptions. The 
need to make such choices introduces a certain bias into the investigation, which can be 
somewhat unsatisfying, especially when hidden or exotic correlations are present. 



1.1. The correlation density matrix 

The correlation density matrix (CDM) [1] has been proposed as an unbiased tool to 
discover the dominant kind of correlations between two separated clusters, given the 
density matrix for their union (obtained by tracing out the rest of the system). For two 
disjoint, separated clusters A and B the CDM is defined to be the density matrix of 
their union minus the direct product of their respective density matrices to get rid of 
trivial correlations, 

pC^pAuB_^A^^B^ (1.1) 
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which is completely unbiased except for the specification of the clusters. If the two 
clusters were not correlated at all, this would imply p^^ = p"^ ® and therefore 
= 0. The CDM encodes all possible correlations between the clusters A and -B, as 
can be seen from the fact that 

tr [p^O^ ® O'^) = tr (p^^^ (6^ ® 6'^)) - tr ((p^^O^) ® (p^6'^ 

= {d^d'^)-{d^){d'^)^Coo>, (1.2) 

where O"^ and O'^ are operators acting on clusters A and B, respectively. 

1.2. Lessons from Luttinger liquid theory 

To extract useful information from the CDM, it will be helpful to develop some intuition 
for its general structure. To this end, let us recall some fundamental facts from one- 
dimensional critical fermion systems. They are described by Luttinger liquid theory, in 
which one of the key parameters is the Fermi wave vector k-p . The asymptotic behavior of 
any kind of correlation or Green's function is typically an oscillation inside a power-law 
envelope, 

C (r) ~ cos {mkpT + (p) /r"' , (1-3) 

for some exponent 7, where m is some integer. For the particular model to be used in 
this study, a nontrivial mapping is known to a free fermion chain [2], a special case of 
Luttinger liquid. 

Renormalization group theory [U] quite generally implies the existence of scaling 
operators in any critical system such as a Luttinger liquid. They are eigenvectors of 
the renormalization transformation and consequently their correlations are purely of a 



form like (1.3) for all r, not just asymptotically. The scaling operators usually have 
complicated forms. The correlation of a simple operator (e.g. fermion density n{x) at 
position X along a chain) has overlap with various scaling operators, and correspondingly 
the correlation function of that simple operator is a linear combination of contributions 



like (1.3) from those scaling operators. 

Our aim is to discover the leading scaling operators numerically. The leading scaling 
operator encodes all the local fluctuations that are correlated with faraway parts of the 
system. Intuitively, for a given cluster A, that operator does not depend signiflcantly 
on the exact position of the (distant) cluster B. That is particularly obvious in a one 
dimensional system: any correlation at distances r' > r must be propagated through 
some sort of correlation at r, so we expect the same operators from cluster A to be 
involved in p*" (r), irrespective of the distance r. 

This suggests an ansatz for leading contributions in the CDM: 

iksr 



f{r) = Y,0^''®0^''cs—. (1.4) 



Here O^''^ and O^'** are a pair of (distance-independent) scaling operators acting on 
clusters A and B, respectively, kg is the characteristic wave vector for oscillations in their 
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correlation, and 7^ is the corresponding scaling exponent. When kg ^ 0, the operator 
pairs must themselves come in pairs, labelled, say, by s and s + 1, with fc^+i = —kg, 
Cs+i = c*, and 7s+i = 7s, so that p*" is hermitian. The scahng operators for each cluster 



form an orthonormal set. We expect that only a few terms in the sum in (1.4) capture 
most of the weight. Correspondingly, it may be feasible to truncate the complete basis 
sets O^'* and O^'* to a smaller set of "dominant operators", whose correlators carry 



the dominant correlations of the system. The ansatz (1.4) will guide our steps in the 



ensuing analysis; at the end, we shall check how well it is satisfied by the actual CDMs 



calculated for the model studied in this paper (see section 6.1.2). 

Notice that although a particular correlation function may have nodes, see (1.3), 
for a CDM of the form ( |1.4[ ) the norm, 

iip'^Mf = E^. (1-5) 

s 

is monotonically decaying with r. This expresses the fact that information can only be 
lost with increasing distance, never restored, in a one-dimensional system. 

1.3. Operator basis and f-matrix 

In [1] the operators entering the dominant correlation were found by a kind of 
singular value decomposition (SVD), which was done independently for each separation. 
However, the operators obtained from the SVD will in general be different for different 



separations r. This does not correspond to the form (1.4), where the operators are 
distance-independent and only the coefficients are r-dependent. Therefore, we shall 
explore in this paper a new scheme to decompose the CDMs for all separations in 
concert, so as to obtain a small set of scaling operators characterizing the dominant 
correlations at any (sufficiently large) separation. We decompose p'" in the form 

p^(r) = ^ (^/'^''^'(r)6^''^®6^''^' I , (1.6) 
where the Si represent the symmetry-sectors of the discrete, Abelian symmetries of 



the Hamiltonian (see section 3.3). The subscript of the brackets indicates that the 
decomposition within the brackets is done for each symmetry-sector individually. This 
decomposition is possible for any two complete, r-independent operator sets O^'^ 
and O^'^' acting on the part of the Hilbert space of clusters A and B, respectively, 
which correspond to the symmetry sector Si. The goal is to find two operator sets 
O^'^ and O^'^' such that these operator sets may be truncated to a small number of 
operators each, while still bearing the dominant correlations of the system. The distance 
dependence of the CDM is then only contained in the matrix f^'^' (r). Then, all analysis 
concerning the distance-dependence of correlations can be done in terms of this f-matrix. 
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1.4- Ground state calculation with DMRG 

The CDM in [T] was calculated using the full ground state obtained from exact 
diagonalization. This limits the system size, so that the method was appropriate mainly 
in cases of rapidly decaying, or non-decaying correlations - not for critical or slowly 
decaying ones. In the present work, we use the density matrix renormalization group 
(DMRG) |3] (see the excellent review by U. Schollwock [1]) to compute the ground 
state for a ladder system which is known to have algebraic correlations [2]. We use the 
matrix product state (MPS) formulation of DMRG [5] in which an efficient variational 
procedure is used to obtain the ground state. 

1.5. Structure of the paper 

The structure of the main body of the paper is as follows: in section [2] we introduce the 
model to be considered for explicit calculations. In section |3] we show how the CDM is 
defined, how to calculate it, and explain how a first overview of the relative strengths 
of various types of correlations can be obtained. In section |4] we show how to analyze 
the CDM and its distance dependence. Sections [5] to [7] present our numerical results, 
and section [8] our conclusions. In an extended appendix we offer a tutorial introduction 
to the MPS formulation of DMRG, and also explain how it can be used to efficiently 
calculate the CDM. 



To be concrete in the following analysis of the CDM, we begin by introducing the model 
for which we did our numerical calculations. This model contains rich physics and its 
treatment below can readily be generalized to other models. 

2.1. Definition of the model 

We analyze the CDM for a class of spinless extended Hubbard models for fermions, which 
was intensely studied by Cheong and Henley [2]. They computed correlation functions 
up to separations of about r = 20, using nontrivial mappings to free fermions and 
hardcore bosons. The correlation functions are calculated with an intervening-particle 
expansion |2], which expresses the correlation functions in terms of one-dimensional 
Fermi-sea expectation values (an evaluation of the CDM for that model has also been 
done by Cheong and Henley P, using exact diagonalization, but the system sizes are 
too short to be conclusive). For spinless fermions on a two- leg ladder with length A^, 
we use the following Hamiltonian: 



2. Model 



2 N-l 
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a=l x=l 



x=l 
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Figure 1. Ladder model with the terms of the Hamiltonian m (2.11. Fermions are 



depicted by black circles and empty lattice positions by white circles. The ordering used 
for our Jordan- Wigner transformation of fcrmionic creation and annihilation operators 
is depicted by the red line. 



2 Af-l N 

+ ^ ^ ^ ^a,x^a,x+l + ^ , (2.1) 

a=l x=l x=l 

where Ca^x destroys a spinless fermion on leg a and rung x, and ha,x = cl^x^a,x is the 
corresponding number operator. Effectively, the model corresponds to a one-dimensional 
pseudo-spin chain, where the a = 1 leg is denoted by spin | and the a = 2 leg by spin 
|. Hence, in the following sections which generally apply to quantum chain models we 
will treat this model as a quantum chain consisting of sites and return to view the 
system as a ladder model in the sections where we discuss our results. 

We will focus on infinite nearest-neighbour repulsion V oo, which we treat 
differently along the legs and the rungs in our numerical calculations. In the pseudo- 
spin description we can enforce the nearest-neighbour exclusion along rungs by removing 
double occupancy from the local Hilbert space of the pseudo-spin sites. The nearest- 
neighbour exclusion along the legs cannot be implemented so easily and we mimic 
— > CX3 by a value of V which is much larger than all the other energies in the 
Hamiltonian (typically V/t\^ = 10^). 

For fermionic systems, the fermionic sign due to the anti-commutation relations 
of the fermionic creation- and annihilation-operators needs to be taken into account. 
Specifically, we have to choose an order in which we pick the Fock basis, where we 
have to keep in mind that this choice produces a so called Jordan-Wigner-string of the 
form '^x'~=x+i e*'^"^" when evaluating correlators {cxcl,) at distance r = \x — x'\. In the 
present system it is convenient to choose this order such that the operators of the two 
sites of a rung are succeeding each other (see figure [T]), as this choice yields the shortest 
Jordan- Wigner strings. 



2.2. Expectations for simple limiting cases 

Setting t|| = 1 as a reference scale, we are left with two parameters in the Hamiltonian: 
the rung hopping tj_ and the correlated hopping t^. The physics of the system is governed 
by the competition of t± to localize the fermions on the rungs and tc to pair the fermions. 
There are three limiting cases which have been studied in detail by Cheong and Henley 

(i) The paired limit, t^ ^ (we used tcA|| = 10^ and t_L = for our calculations). 
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In this limit the fermions form tight pairs which behave similar to hardcore bosons 
[2]. For two given rungs x and x + 1, there are two possibilities to create a pair of 
fermions, due to infinite nearest-neighbour repulsion: c|^c|^_|_^ and c|2,c|^._^]^. It has 
been shown in |2] that, based on these two bound pairs, one may classify the bound 
pairs in two flavours along the ladder and that the ground state has only one definite 
flavour, causing a twofold symmetry breaking in the ground state. This symmetry 
breaking introduces complications that will be addressed below. The dominant 
correlations are expected to be charge-density correlations at short distances and 
two-particle at long distances. These charge-density and two-particle correlations 
decay as power laws, oscillating with k = 2kF, where the Fermi wavelength kp is 
related to the filling as fcp = 2i/ [2]. In this system, the one-particle correlations are 
suppressed and are expected to decay exponentially, as a nonzero expectation value 
depends on a local fluctuation completely filling the rungs between the clusters (as 



elaborated in section 6.2). 

(ii) The two-leg limit, tj_ -C ty, t^. = 0. In this limit the two legs are decoupled with 
respect to hopping, but still the infinite nearest-neighbour repulsion introduces 
correlations between the two legs. At large distances, power-law charge- density 
correlations dominate, while two-particle correlations show much faster power-law 
decay and one-particle correlations decay exponentially. 

(iii) The rung-fermion limit, t± ^ t\\, tc = 0. In this limit the particles are delocalized 
along the rungs. For fillings smaller than quarter-filling, charge-density , one- 
particle and two-particle correlations all decay as power laws where charge- density 
correlations dominate at large distances. 

Our analysis in this paper is limited to the case (i), where DMRG also showed best 
performance. 



2.3. Smooth boundary conditions 

For a ladder of length (treated as a pseudo-spin chain) , we have attempted to reduce 
effects from the boundaries by implementing smooth boundary conditions, adopting a 
strategy proposed in [7j for a spin chain to our present fermionic system. (Alternatively, 
it is possible to use periodic boundary conditions [5]. However, this leads to some 
difficulties, since it is not possible to work with orthonormal basis sets describing 
the left or right part of the chain with respect to a given site.) Smooth boundary 
conditions are open boundary conditions together with an artificial decay of all terms 
of the Hamiltonian over the last M rungs at each end of the chain. We shall calculate 
expectation values only of operators located in the central part of the system (sites x, 
with M < X < N - M), thus the system's effective length is A^' = A^ - 2M. 

For both smooth and open boundary conditions the average site filling strongly 
decreases near the boundaries. To determine the average filling z/, which influences the 
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system's correlations in an important manner, we thus use only the central A^' sites: 



N-M 



y= E ((^T.) + (^ix))/(2iV')- (2.2) 
Due to the infinite nearest neighbour repulsion, this implies that v G [0,0.5]. 

3. Calculation of the CDM 

Throughout the paper we will use the Frobenius inner product and norm for any matrices 
Mij and M[- of matching dimension, 

(M, M') = J2 ^h^ij = [M^M') (3.1) 
||M|| =(M,M)^/^ (3.2) 
3.1. Definition of the CDM 

We take two disjoint, separated clusters A and B of equal size from a one-dimensional 
quantum chain, i.e. two sets of adjacent sites x^, . . . and a;f , . . . , where n is the 
size of the clusters and all the indices x are distinct from each other. The local Hilbert 
spaces of clusters A and B with dimension c?" are described in terms of sets of basis 
states I a) and |/?), which are product states of the local states of each site in the cluster. 



The CDM of the two clusters, defined by (1.1), can be expanded in this basis as 

p"" = P%.,p,\a)\P){a'\m . (3.3) 
For processing the CDM we fuse the two indices of each cluster 

P%^PUw)\^){^'m{(^'\ (3.4) 

with a = {aa') and (3 = {13(3'), and denote the reshaped object itself by an extra 
tilde. This corresponds to a partial transpose of the CDM (note that p*" is no longer 
a symmetric tensor). For the CDM expressed in the indices a and /5, we may use the 



Frobenius inner product (3.1) and norm (3.2). 



To study the distance dependence of the correlations, we vary the position of the 
clusters A and -B, resulting in a position-dependent CDM p'^ (a;^, xf ). If the system is 
translationally invariant, this object depends only on the distance r = — xf | (the 
minimal distance for two adjacent clusters is equal to the cluster size n). For a finite 
system, though, p^ will also depend on | [x^ + xf ), at best weakly if the system is long. 
Strategies for minimizing the dependence on | [x^ + xf ) by taking suitable averages 
will be discussed in section 13.41 



3.2. DMRG- calculation of the CDM 



The fact that the Hamiltonian in (2.1) is a one-dimensional pseudo-spin chain allows 
us to calculate ground state properties with the density matrix renormalization group 
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(DMRG) [3]. Using the variational matrix product state formulation of that method (see 
appendix for a detailed description), we calculated the ground state of the Hamiltonian 
in (2.1) for several values of t± and tc- The framework of MPS also allows the CDM 
to be calculated efficiently (see section A. 2. 7 for details). Limiting ourselves to the case 
t± = in this paper, we have calculated the CDM derived from the ground state for 
distances up to 40 rungs, which is significantly larger than in previous approaches. 



3.3. Symmetry sectors 

All the symmetries of the Hamiltonian are refiected in the CDM, making the CDM 
block-diagonal, where each block can be labeled uniquely by a set of quantum numbers 
that are conserved by the Hamiltonian. This means for Abelian symmetries (which 
are the only ones we are considering in practice), that the CDM in the original form 
Pa/3,«'/3' fulfills Qa + Q /3 = Q a' + Q 13' ^ where Qa correspouds to the quantum numbers 

then implies AQ^ = — AQ^ 

of the CDM involving AQq, (AQ^) there has to be a block involving —AQ^ {—AQ^), 
respectively. Therefore, it is convenient to sort the various parts of the CDM in terms of 
their change in quantum numbers AQ = \AQa\ = \AQ^\ and to analyze each symmetry 
sector individually. 

To obtain a general classification of the CDM we sort the various contributions 
of the CDM according to the conserved quantum number(s) Q. In the case of the 



of state \a), etc. The rearrangement of the CDM into p?^ then implies AQ^ = —AQj. 
with AQa = Qa — Qa' and AQa = Qp — Qpi. Since p^'^ is hermitian, for every block 



Hamiltonian in (2.1), we consider particle conservation {Q = Ntot) which breaks the 
CDM into blocks with well-defined particle transfer AA^ = lAA"^! = |AA^g| between 
clusters A and B. The following r.m.s. net correlations then is a measure for the 
correlations with transfer of AA^ particles between A and B (with AA^ = 0, 1, 2): 

<N{r)= E (3-5) 

where X]AAf=o "^atv (^) — llp'"(^)lP- Here the notation a = {aa') G San indicates 
that only pairs of states (aa') are considered which differ by AA^ in particle number 
(similarly for (3 = {(313') G San)- In the following we will call correlations involving 
AA^ = 0,1,2 particles charge- density correlations (CD), one-particle correlations (IP), 
and two-particle correlations (2P), respectively. The following analysis is done for each 
symmetry sector individually. Depending on the decay of the r.m.s. net correlations 



(3.5), some symmetry sectors may become irrelevant with increasing distance. 



3.4. "Restoration" of numerically broken symmetries 

Although we have tried to minimize the effect of boundaries, our numerical methods for 
calculating the ground state and CDM do not produce strictly translationally invariant 
results. (In contrast, analyses based on exact diagonalization start from a ground state 
wavefunction in which the symmetry (in a finite system) is restored, even if there is 
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a symmetry breaking in the thermodynamic hmit.) Therefore, we construct the CDM 
p*-^ (r) for a given distance r from an average over several CDMs p'" {x, x') with constant 
r = |x — x'l, where x and x' give the position of the first site of clusters A and B, 
respectively. 

Moreover, if the exact ground state is degenerate under a discrete symmetry, we 
expect that DMRG breaks this symmetry unless it is implemented explicitly in the 
code. As mentioned in section 2.2 for the specific models of this paper we expect a 
discrete symmetry under interchange of legs for some parameter regimes. Since we 
did not implement this symmetry explicitly in our code, we also average the CDM by 
interchanging the legs of the ladder. Thus, all the data analysis presented in subsequent 
sections will be based on using the following "symmetry-restored" form of the CDM, 

p^ (r) = 1 J2 ('P"" + 'P"" ^')) , (3.6) 

xx' ,\x—x'\=r 

where p"^ is obtained from p^ by interchanging the legs of the ladder, and A/" is some 
normalization factor. 

One might argue that it is not sufficient to average over the broken symmetry w.r.t. 
leg-interchange on the level of the density matrix, but that instead the symmetry should 
be restored on the level of the ground state wave function. Specifically, for a ground 
state \ipi) (however it is calculated) which breaks this symmetry, we could restore the 
symmetry in the following way, 

|^+> = ^(lV^i) + l^2)) , (3.7) 

where \1IJ2) = S {ipi) and S describes the action of interchanging the legs. This would 
lead to a total density matrix 

li^-") {i^-^l = I (IV^i) (^il + \H iH + li^i) iH + 1^2) (V^il) . (3.8) 



Now, for two clusters A and B, the first two terms on the r.h.s. yield the CDM of (3.6), 
while the last two terms turn out to be negligible when traced out over all sites except 
for the two local clusters A and B. This follows from lipi) and \1IJ2) being orthogonal, 
hence tT{\ijji) (V'2|) = {'4'2\'ipi) = 0, implying that for a long chain with local clusters A 
and B, the reduced density matrix p^-^'^^ ^ tr^,^^ b(|'?/'i) (V'2|) will be very close to zero 
due to the orthogonality of the wave functions on the sites outside of clusters A and B. 
Consequently, it is sufficient to retain only the first two terms of (3.8), i.e. to restore 
the broken symmetry on the level of the density matrices only, as done in (3.6). 



4. Finding a distance-independent operator basis 

The goal of this section is to extract a (likely) small set of operators from the CDM, 
which will describe the dominant correlations in the system as a function of distance. 
We will assume in this section that the CDM does not include any broken symmetries 
as indicated in section 13.41 
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4.I. Need for operator bases for clusters A and B 



As already mentioned, the CDM (obtained from (3.6)) may be investigated by applying 
a singular value decomposition (SVD) for each distance individually P: 

~p%=Y.'''or®of\ (4.1) 

s 

or, in operator notation: 

s 

where O^''^ and O^'^ act on clusters A and -B, respectively. Here the singular values 
are strictly positive real numbers. By construction, O"^'* and O^''^ form orthonormal sets 
in their corresponding Hilbert spaces, i.e. O^'"* = O^'^, and O^'^ = O^^, form a complete 
set in the operator space of clusters A and B, respectively, using the inner product as 



in (3.1). The set includes operators with Ws = 0, such as the identity operator, since 



these will be produced by the SVD. The SVD (4.2) yields for each specific distance r a 
set of operators O"^'* (r) and O^''^ (r) acting on clusters A and B, respectively. 

However, the dominant operators so obtained, i.e. the ones with large weight from 
the SVD of p*" (r), are likely not the same as each other for different distances and hence 
not convenient for characterizing the "dominant correlations" of the system. What is 
needed, evidently, is a strategy for reducing the numerous sets of operators O^'^ (r) and 
qb,s "basis sets of operators" for clusters A and B, respectively, say O^''^ and 

O^'^, which are r-independent and whose correlators yield the dominant correlations in 



the system in the spirit of (1.4). (For a translationally invariant system the two sets 



have to be equal for both clusters A and B, but we will treat them independently in the 



analysis.) Following the ansatz (1.4) from the Luttinger liquid theory, these operators 
ought to be distance- independent, carrying common correlation content for all distances. 
Thus we seek an expansion of (r) of the form (1.6), in which only the coefficients, 
not the operators, are r-dependent. 



4-2. Construction of operator bases 

We have explored a number of different strategies for extracting operators from the 
CDM which carry common information for all distances. We will discuss in detail only 
one of these, which is rather simple to formulate and reliably yields operator sets with 
the desired properties. (Several other strategies yielded equivalent results, but in a 
somewhat more cumbersome fashion.) 

The simplest possible strategy one may try is to average over all the CDMs at 
different distances and to singular- value decompose the resulting crude "average CDM" . 
However, since the elements for the CDM are expected to be oscillating functions of r, 
such a crude average can cancel out important contributions of the CDM. Thus we 
need a procedure that avoids such possible cancellations. To this end, we construct the 
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following operators, bilinear in the CDM: 

i^^(r) = tiB{fHr)f{r))/\\f\\' (4.3a) 
k^{r) ^ tTA{f{r)fHr))/\\fr, (4.36) 
with matrix elements 

/3 

We normalize by ||p'^(r)||^ in order to treat the operator correlations of {r) for 
different distances on an equal footing. Note that the eigenvalue decomposition on the 
hermitian matrices K"^ (r) and (r) (in short K-matrices) yields the same operators 
O"^ (r) and (r) as the SVD of p'" (r), with eigenvalues being equal to singular values 
squared, up to the additional normalization factor ||p*" (r) |p. (Reason: for a matrix of 
the form M = usv^ we have MM'' = us'^u' and M'^M = f s^f ^.) 

The object (for X = A, B) is positive-definite and according to ansatz (1.4), it 
is expected to have the form 

K''ir)=Af,^J2\^d-d-K (4.5) 



In particular, it no longer contains any oscillating parts (in contrast to (1.4)), and hence 
is suitable for being averaged over r. 

Summing up the i^^^-matrices over a range R of distances (r G R, where R 
will be specified below) gives a mean ^"^-matrix for cluster X (= A,B), namely 
^x,R ^ ^^g^ (r). We do not divide the latter expression by the number of terms in 
the sum (as would be required for a proper mean), as at this stage we are only interested 
in the operator eigendecomposition, 

with the operators normalized such that ||0"^''^''^|| = 1. The operator set O^'^''^ gives 
an orthonormal, r-independent basis for cluster X. In practice, however, many of the 
^■R.M (which turn out to be the same for X = A or B) will be very small. Thus, it will 
be sufficient to work with a truncated set of these operators having significant weight. 

To explore the extent to which depends on the summation range, we shall 
study several such ranges: i?aii includes all distances, i?short short distances (first third of 
distances analyzed), i?int intermediate distances (second third) and -Riong long distances 
(last third). The resulting (truncated) sets of operators can be compared via their 
mutual overlap matrix O^^' = tT{0^'^''^0^''^''^'), or more simply, by the single number 
qBR' _ ^^^,(0^^ )^, which may be interpreted as the dimension of the common 
subspace of the two operator sets. The value of O^^' ranges from to dim{0^''^'^). 
By comparing O^^' for the different distance ranges, additional clues can be obtained 
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about how the relative weight of correlations evolves from short to long distances. (Such 
a comparison is carried out in table [l] below.) 

4-3. Definition of f- Matrix 

Once a convenient basis of operators O^'^ and O^'^ has been found, the correlation 



density matrix can be expanded in terms of this basis as in (1.6) 



fin' 



with matrix elements 

For complete operator spaces O^'^ and O^'^' , by definition, the set of amplitudes squared 
sum up to the norm of the CDM: 

^|/My(^)|2=||pC(^)||2_ (49) 

However, as alluded to above, we expect that the dominant correlators can be expressed 
in terms of a truncated set of dominant operators. If the sum on the left hand side of 



(4.9) is restricted to this truncated set, its deviation from the right hand side gives an 
estimate of how well p'" is represented by the truncated set of operators. It will turn 
out that only a handful of dominant operators (typically 4 or 6) are needed, implying 
very significant simplifications in the analysis. Thus, the data analysis will be done in 
terms of the matrices f^'^' (r) (in short "f-matrix" ) for this truncated set of dominant 
operators. 

4.4- Fourier- analysis and decay of f-matrix 



According to the expectations expressed in (1.4), the elements of the f-matrix are 
expected to be products of oscillating and decaying functions of r. The corresponding 
dominant wave vectors can be identified via Fourier transform on each element of the 
f-matrix. For an oscillating function times a monotonically decaying envelope, the peaks 
of the Fourier spectrum of the oscillating function will be broadened by the presence 
of the envelope. To minimize this unwanted broadening, we introduce a rescaled f- 
matrix (denoted by a tilde), f^^'^^' [r) = u{r) f^'^' (r), where the positive weighting- 
function u{r) is chosen such that all values of \f^'^' (r) | are of the same order, and 
Fourier decompose the rescaled /-matrix as f'^''^' (k) = e"*'^''/^'^' (r). Its norm 
11/ (A;) IP = J2fifi' I/'''''' (^) plotted as a function of k, will contain distinct peaks 
that indicate which wave vectors characterize the dominant correlations. Subsequently, 
the elements of the f-matrix, can be fitted to the forms 
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where A^ji^i are complex amplitudes, fj{r) describes the decay with distance (e.g. 
fj{r) = r~'^i or e~'"/'~J for power-law or exponential decay, respectively), and kj is a 
set of dominant wave vectors. The latter appear pairwise in combinations {+k; —k), 
since f'^''^' G R, which implies ^|^'^/ = fo^^ = The results of such a fit for 
each pair of dominant operators O^''^ and O^'^', is the final outcome of our analysis, 
since it contains the information needed to check the applicability of ansatz (1.3). 



5. Numerical results: general remarks 

In this section, we illustrate the analysis proposed above for the model introduced in 
section|2} We will focus on the limiting case of large tc, which we expect to have the most 
complex behavior among all three limiting cases introduced in [1] and [2]. After some 
preliminary analysis, we will discuss in section |6] each of the three symmetry sectors 
(CD, IP, and 2P) characterized by the operators' fermion number, and in section [T] 
compare our results to those found by [2] using a different method. 



5.1. Specification of the clusters A and B 

For the following analysis it is convenient to take the size of the clusters A and B to be 
two rungs, because clusters of at least that size allow for up to two particles in one cluster 
(due to infinite nearest-neighbour repulsion). Thus, correlations involving AA^ = 0, 1, 2 
are possible, i.e CD, IP, and 2P correlations, respectively. Note that larger clusters 
can be studied, but would significantly increase numerical costs. Taking into account 
the infinite nearest-neighbour repulsion, clusters of size two have a seven-dimensional 
Hilbert space spanned by the kets |00), |0 T), |0 i), |T 0), |i 0), |Ti), lit), where the 
first (second) entry corresponds to the first (second) rung, represents an empty rung 
and I and | a fermion on the upper and lower leg in pseudo-spin notation (recall that 
we are dealing with spinless fermions). The space of operators acting on a cluster has 
dimension 7^ = 49, where the subspaces for AA^ = 0, 1 or 2 have dimensions 21, 24 and 
4, respectively, as depicted schematically in figure |2] 



5.2. Average site occupation 

As a first check of the infiuence of the boundaries, we investigate the average site 
occupation on the ladder. It is expected to be uniform in a translationally invariant 
system. However, there are two ways in which our calculation breaks translational 
symmetry, which cause residual oscillations in the density of particles along the ladders. 
Firstly, there is the spontaneous breaking of the pair fiavor symmetry described in 



section 2.2 In the ground state produced by DMRG, all pairs have the same fiavor, so 
only one of the two sublattices actually has any fermions on it. Thus a strong alternation 
in the density is observed between one leg for even rungs and the other leg for odd rungs; 



this can be taken care of by the symmetrization with respect to legs (as in (3.6)). 



Numerical results 



16 



/ 



|00)|()T>|0i)|T0)|i0>|n>|iT> 



(00 1 
(T0| 
(io| 
(UTI 
{Oil 
UTI 
(Til 
















I 
















h 
















— i 








































2- 










— ( 


)— 



Figure 2. The symmetry sectors of an operator acting on a cluster of two rungs in 
the basis |00), |0 t), |0 j), |t 0), |i 0), |ti), lit) in pseudo-spin notation. 



Secondly, translational symmetry is broken due to finite size in the DMRG 
calculation. This induces oscillations in the average occupation as a function of x (see 
figure [3]), whose period is clearly dependent on the filling. In fact, their period is 2k-p, so 
they may be interpreted as Friedel-like oscillations caused by the boundaries. Although 
the amplitude of density oscillation appears rather flat in the central portion of the 
system, it does have a minimum there; so we expect that the amplitude in the center of 
the system would vanish in a sufficiently large system. 

Although the intent of the smooth boundary conditions is to minimize effects such 
as these oscillations, in fact, their amplitude appeared to be of about the same strength 
independent of whether we used smooth or plain open boundary conditions. We suspect, 
however, that the amplitude could be reduced by further careful optimization (not 
attempted here) of the parameters of the smooth boundary conditions. 



5.3. r.m.s. net correlations w {r) 

The next basic step is to identify the leading correlations in terms of the r.m.s. net 



correlations watv defined in (|3.5|). These reveal which sectors of correlations dominate 

net correlations 



at large distances. The results (see figure |4|) show that the r.m.s. 
decay exponentially in the IP sector, whereas they decay algebraically in both the CD 
and 2P sectors, consistent with |2|- The latter two correlations are comparable in size 
over a significant range of distances, but for the fillings we investigated, 2P correlations 
ultimately dominate over CD correlations at the largest distances. Both the CD and 
2P r.m.s. net correlations can be fitted to power laws, with the exponent dependent on 
the filling. The r.m.s. net correlations in each sector are monotonic and only weakly 
modulated, even though the dominant correlation functions and the dominant parts of 
the CDM itself are oscillating (as will be discussed in more detail in section 6.1, see, 
e.g., figure [?]). This implies that the correlations in each sector can be represented 
by a linear combination of correlation functions (associated with different operators) 
which oscillate out of phase, in such a way that in the sum of their squared moduli the 
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Figure 3. The average occupation along the legs of the ladder for a filling of = 0.248 
(panels a,c) and a filling of = 0.286 (panels b,d). Panels (a) and (b) show the average 
occupation nij on the upper (red) and hi on the lower (green) leg, with every second 
value being zero. The end regions « = 1 ... 20 and i = 81 . . . 100 were skipped in the 
figures and also in the analysis as these are affected by the smooth open boundary 
condition. The leg symmetrized occupation h = ^(n| + hi) (blue, same for upper and 
lower leg) eliminates this strong even odd alternation but still shows small modulations. 
This can be seen in detail in the Fourier transform of the symmetrized occupation in 
panels (c) and (d). There is a clear peak at A; = ±2A;f (dashed vertical lines). 



oscillations more or less average out, resulting in an essentially monotonic decay with 
r, as expected according to (1.5). 



We will next apply the analysis proposed in section 4^ to the respective symmetry 
sectors (which will provide more exact fits of the exponents of the power-law decays). 



The analysis in any sector consists of two stages. First, following section we try to 
find an optimal truncated basis which describes best the dominant correlations. Second, 



we examine the f-matrix of section 4.3 (i.e. represent the CDM in the truncated basis) 
to see the nature of its r dependence, and to fit this to an appropriate form, following 
section I 



6. Numerical results: symmetry sectors 
6.1. Charge- density correlations 

6.1.1. Operator basis First we calculated the mean K-matrices K"^'^ and K^'^ from 



defined in (|4.3a[) and (|4.36[), and obtained operator sets from their eigenvalue 
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Figure 4. The r.m.s. net correlations of (3.51, plotted as a function of distance for (a) 
a filling oi V — 0.248 and (b) a filling oi v ~ 0.286. The symmetry sectors are AA^ = 
(blue, no particle transfer, CD), AiV = 1 (green, transfer of one particle, IP) and 
AA^ = 2 (red, transfer of two particles, 2P). We see that CD and 2P correlations decay 
as power-laws {r~'^ , blue and red solid lines) with small residual oscillations at fc = 2fcFj 
while the IP correlations show exponential decay [e^^/^^ , see semi-logarithmic plot in 
the inset). The value ri ~ 0.5 for both fillings is reasonable as we would expect a value 
of the order of one, which is the size of the bound pairs. 



decomposition, using various distance ranges. 

In order to decide how many operators to include in the truncated basis, we used 
the diagnostic described in section 4^ In presenting the results, we limit ourselves 
to cluster A as the results for cluster B are completely analogous. The operator set 
QA,R^ii,ti corresponding to the full range of distances i?aii (specified in section section 4.2) 
is used as a reference set to be compared with the operator sets obtained from -Rghort, 
i?int and -Riong- The results are given in table [ij We see that, for intermediate or long 
distances, the effective dimension (O'^''"^'"* and 0^='"'^'°"s) of the common operator space 
shared between the operator set O"^'-^''"''^ and the operator sets O^^'-^^"*''^ and 0^'^i°ns''^, 
respectively, saturates at six even if a larger operator space is allowed. Similarly, also 
the short-distance operator set agrees best with the other three operator sets 

at dimension six: a further increase of the number of operators, however, adds only 
operators in the short range sector of the CDM. Hence we truncate to a six-dimensional 
operator basis. Within this reduced operator space, all dominant correlations are well- 
captured, as can be seen from the relative weights of table [Tj For the resulting truncated 
basis set equation (4.9) holds up to a relative deviation of the order of O (10~^). 

Investigating the six-dimensional set of operators in more detail reveals that they 
can be classified with respect to their symmetry with respect to interchanging the legs 
of the ladder, i.e. they obey 5*0^'^='"''^ = iO^'^aii,/^^ with S describing the action of 
interchanging legs. The set breaks into two subsets of three operators each, which have 
positive or negative parity with respect to S, respectively. It turns out that all six 
operators are linear combinations of operators having matrix elements on the diagonal 
only, in the representation of figure [2j Moreover, together with the unit matrix they span 
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Table 1. Comparison of the operator sets on cluster A for a filling of = 0.286. (The 
results for v = 0.248 are similar, with only minor difl^erences.) The first and second 
column of the table give the number of operators kept and the corresponding smallest 
singular value of the set of operators O"^^^""-'^ obtained from the full range of distances 
i?aii- The other three columns show O^'^"^'^^"* , O^"^"^^"* and O^^iifliong fQj. ^i^q given 
number of operators. 



number of jy^Ra.W-'^ (^^all^short ^^all^int (^-Rall -Rlong 

operators (short) (intermediate) (long) 



1 


1 


1 


0.99 


1 


2 


0.784122 


1.99 


2 


2 


3 


0.579242 


2.99 


3 


3 


4 


0.176043 


3.99 


4 


4 


5 


0.011250 


5 


5 


4.99 


6 


0.003040 


6 


6 


5.99 


7 


0.000004 


7 


6 


6 


8 


0.000001 


8 


6 


6 


9 


0.000001 


9 


6 


6 


10 


0.000001 


10 


6 


6 



the full space of diagonal operators (therefore the dimension of 6 = 7 — 1). Explicitly, 
the symmetric operators are given by 

6^ = ^ {-fiQ^xfi^^x+i - n^,xno,x+i + 2h^^xh[,x+i + leg symmetrized) (6.1a) 



0=2 {no,xn^,x+i - n'i,xno,x+i + leg symmetrized) (6.16) 
= ^ [-Gho^xfio^x+i + {fio^xfi^^x+i + n'^,xno,x+i + n^,xni,x+i + leg symmetrized)] (6.1c) 

and the antisymmetric operators by 

= Ts'^o.x {n^,x+i - fii^x+i) (6.2a) 
= ^ {n^,x - ni,x) no,x+i (6.26) 
= ^ {h'^,xni,x+i - fii^xfi^^x+i) (6.2c) 

where no = (1 — rif — fii). We use this operator basis for both cluster A and cluster B. 
If we calculate the f-matrix (4.7) based on these operators we see that it breaks into two 
blocks corresponding to their symmetry with respect to leg interchange. 



6.1.2. f-matrix elements: oscillations and decay We now turn to extracting the 
distance-dependence of the dominant correlation in this symmetry sector, which is 
now visualizable since we drastically reduced the operator space to six dimensions. 
All relevant information is contained in the f-matrix and its Fourier transform. The 
first step is to identify the oscillation wave vector(s) k to be used as initial guesses in 
the fit. A general method is to plot the Fourier spectrum ||/(fc)|| of the rescaled f- 
matrix (figure |5]). When using a logarithmic scale for the vertical axis, even sub-leading 
contributions show up clearly. We find that the spectra belonging to the symmetric 
and anti-symmetric operators are shifted against each other by vr. This relative phase 
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Figure 5. Fourier transform of the rescaled f-matrix / for CD correlations based 
on operators chosen from a reduced six-dimensional operator space, for a filling of 
(a) v = 0.248 and (b) v = 0.286. We obtain these Fourier spectra from the 
rescaled f-matrix /^'^ (r) = r'^ /^■'^ (r), with 7" extracted from a power-law fit 
on 1/'^''^ (r) |. The Fourier spectrum breaks up into a contribution coming from the 
operators symmetric or antisymmetric under leg-interchange, labelled (blue) and 
/~ (red), respectively. The spectrum of /+ shows strong peaks at fc = ±2fcF (dashed 
lines) and a smaller peak at fc = with kp/ir = v. The spectrum of /", having peaks 
at fc = ±2fcF -I- TT (dashed lines) and fc = tt, is shifted w.r.t. /+ by tt. For a filling close 
to i the dominant peaks of at fc = ±2fcF and fc — ±2fcF + t^- are nearly at the 
same position. 



shift implies a trivial additional distance dependence of e™"^ of f~{r) with respect to 
/"•"(r), reflecting the different parity under leg interchange of the two operator sets. We 
have found it convenient to undo this shift by redefining f~{r), the part of the f-matrix 
belonging to the anti-symmetric operators, to e*'^'"/" (r). The resulting combined Fourier 
spectrum for /"*" and e*'^''/" has strong peaks at = 2kY and a smaller peak at fc = 0, 
in agreement with the result from [2J. 

Based on the Fourier spectrum, we rewrite the fitting form (4.10) as 

r^'^' (r) = A^^,r-^ cos (fcr + 0^^,) + B^^,r-^' , (6.3) 

with real numbers A^^i > and -B^^^', where we expect 7' > 7, due to the relative 
sharpness of the peaks in the Fourier spectrum. The non-linear fitting over the full 
range of distances is done in several steps to also include the decaying part at long 
distances on an equal footing. First, the data is rescaled by r^^" , where we obtained 7" 
from a simple power-law fit, in order to be able to fit the oscillations for all distances 
with comparable accuracy. Then we fit the rescaled data to (6.3), where initially we use 
the information from the Fourier spectrum in keeping fc fixed to fc = 2fcF, but finally 
also release the constraint on fc. This procedure showed best results, with relative error 
bounds up to 2%. The uncertainties are largest for the second term in (6.3) as it acts 
mainly on short distances, having 7' > 7. 

The results of this fitting procedure are depicted in figure |6| for all 18 nonzero 
elements of the f-matrix. We see that the leading power-law exponents deviate from the 
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Figure 6. The results of the fit in (6.3) to the 18 independent elements /'^■'^ of the 
f- matrix, labeled along the horizontal axis by the index pair fi/i' , for AA'^ = at filling 
h' — 0.286. The results for = 0.248 are qualitatively the same. Panel (a) shows 
7, panel (b) ^/i^', panel (c) 0^^' and panel (d) the error of the fitting e, defined by 
= ^^{f^^ {r) — f^^{r))'^/r~'^^ , where r~'^ is the power-law we used to rescale the 
data before Fourier-transforming. The red, dashed line in the first panel shows the 
power-law exponent obtained from the r.m.s. net correlations, 70 = 1.33. The phase 
(j)fifj_i is defined such that it is in the interval [— 7r,7r]. The matrix elements have been 
grouped according to their relative phases (/)^^' (separated by the black, dashed line), 

= and 0^^' — respectively. 



which clearly indicate cos and sin behaviour for 



The solid red lines in panels (a) and (b) show the exponent 70 and the amplitude A, 



respectively, from the single fit (6.4 1 



fit to the r.m.s. net correlations in (3.5) (compare figurejl]) by about 5%. The /c-vectors 
from the non-hnear fit are close to = 2kp and deviate by less than 1%. The fit to 



the sub-leading second term in (6.3) is not reliable, so we do not show the results for 7' 



here, but note that every fit satisfied 7' > 7. 

Since most of the exponents 7 and amplitudes A^j^^i are of comparable size, we fit 
the f-matrix elements to a single 70 and A (as well as a single 7q and B for the second 
term) for all the f-matrix elements, using the Ansatz: 



/ (r) = Ar 



■70 



V 



cos(A;r) 
— sin(fcr) 
cos(A;r) 



sin(fcr) 
cos(fcr) 
sin(A,T) 





cos(fcr) 
— sin(A;r) 
cos(fcr) 





cos(fcr) cos(A;r) sin(fcr) 
cos(fcr) cos(A;r) sin(fcr) 
~ sin(fcr)— sin(fcr)cos(A;r) 



Br '^u 



(6.4) 



The form of the matrices in the two blocks was obtained by inserting into (6.3) the 
explicit values of the phases 0^^/ determined from the previous fit and summarized in 
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Figure 7. Two entries of the f-matrix for (a) ly ~ 0.248 and (b) v — 0.286 fitted to 



the form in (6.3). The single points (blue circles and squares) are data points from 
the f-matrix and the lines (red and green) are the result of the fitting. They evidently 
oscillate with a relative phase of A(j) = 7r/2. As a result, their contribution to the 
r.m.s. net correlations , (|/^'^P + l/^'^P)', shown by the thick orange curve, has only 
small oscillations at large distances. 



figure |6} Fitting to (6.4) gives an error of about 10%, witli largest errors arising for 
the f-matrix elements where A^^i deviates strongly from A (see figure |6]). For the filling 
V = 0.286 we find 70 = 1.26 and A = 0.06. The values of and B are unreliable in 
that the results from several fittings differ by about 30%, but still it holds that 7o > 7o- 



The form of (6.4) allows us to understand why the r.m.s. net correlations displayed 



in figured show some residual oscillations, instead of decaying completely smoothly, as 



anticipated in section 1.2 The reason is that (6.4) contains 10 cos{kr) terms but only 
8 sin(A;r) terms. Although any two such terms oscillate out of phase, as illustrated in 
figure [7], the cancellation of oscillations will thus not be complete. Instead, the r.m.s. 
net correlations contain a factor [8 + 2cos^(fcr)]^ (compare to (3.5)), which produces 
relative oscillations of about 10%, in accord with figure |4j (The fact that the total 
number of cos(/cr) and sin(fcr) terms is not equal is to be expected: the total operator 
Hilbert space per cluster is limited, and its symmetry subspaces might have dimensions 
not a multiple of 4.) 

For each pair of wave vectors ±k in each parity sector, the effective operator basis 
per cluster can be reduced even further, from 3 operators to one conjugate pair of 



operators. This can be seen by rewriting (|6.4|) as follows: 

/+ 



/ (r) = Ar~ 



70 



Akr 







e^'^V- 

with the matrices /+ and /_ defined as 



+ c.c. 



(6.5) 



/ 1 


-^ l\ 


' 2 


(l 1 


-A 


i 


1 i 


1 1 


—i 




-^ I) 




^ z i 


1 / 



(6.6) 
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Table 2. Comparison of the IP operator sets on cluster A for a filling oi u — 0.286, 
using the same conventions as for table [T] 



number of m^aU^t^ ^^R^n-A (^^all^short ^^all^int (^-Rall-Rlong 

operators (short) (intermediate) (long) 



4 


1 


4 


4 


4 


8 


0.297162 


8 


8 


8 


12 


0.014661 


12 


12 


12 


16 


0.000402 


16 


16 


16 


20 


0.000001 


19.97 


19.95 


19.31 



Note that both and /_ are matrices of rank one with eigenvalues |, and 0. The 
eigenvectors with eigenvalue | are ;^(l,i, 1) and -^{1,1, i), respectively. Thus, by 
transforming to an operator basis in which f± is diagonal, one finds that in both the 
even and the odd sector, the dominant correlations are actually carried by only a pair of 
operators, namely -^{O^ + iO^ + 0^) and its hermitian conjugate, and ■^{O'^ + O^ + iO^) 
and its hermitian conjugate, respectively. This result, whose precise form could hardly 
have been anticipated a priori, is a pleasing illustration of the power of a CDM analysis 
to uncover nontrivial correlations. 



6.2. One-particle correlations 

The correlations in the IP sector are exponentially decaying, as already mentioned in 



section 5.3 The reason for this was given in py and is the key to understanding the 
operators and correlations in this sector. In the limit where the fermions are all paired, 
the only possible way to annihilate one at x and create one at x' > x , such that the 
initial and final states are both paired, is that every rung in the interval {x,x') has a 
fermion (necessarily on alternating legs). These fermions can be grouped as pairs in 
two different ways: {x,x + 1), {x + 2,x + 3), . . . , (x' — 2, x' — 1) in the initial state, 
but (x + 1, X + 2), . . . , (x' — 1, x') in the final state. (Notice this requires that x and x' 
have the same parity.) pQ showed that the probability of such a run of filled sites decays 
exponentially with its length. 

Applying the operator analysis in this sector using the eigenvalue decomposition in 



(4.6) gives a series of fourfold degenerate eigenvalues for both clusters, see table |2j for 
cluster A. The table for cluster B is exactly the same. For a specific eigenvalue, also 
the operators for cluster B (residing at rungs (x',x' + 1)) are the same as for cluster A 
(residing at rungs (x, x + 1)), but with mirrored rungs, i.e. an operator acting on rungs 
(x, X + 1) acts in the same fashion on rungs (x' + l,x'). 

Looking more closely, the first four operators annihilate or create a particle on 
rungs X + 1 or x', respectively, thereby breaking or regrouping bound pairs residing on 
(x + l,x + 2) or (x' — l,x'), respectively. The second set of four operators annihilates 
or creates a particle on rungs x or x' + 1, thereby breaking or regrouping bound pairs 
residing on rungs (x,x + 1) or (x',x' + 1). For a given odd separation x' — x, the 
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Figure 8. Three configurations of bound pairs contributing to IP correlations for a 
distance (a,) r — 2 and (b),(c) r — 3. Clusters A and B are depicted by the green 
and red squares, respectively. Fermions are depicted by black circles, empty lattice 
positions by white circles and the position where a fermion will be created is depicted 
by concentric circles. The crosses show the center of mass of the bound pairs. In 
configuration (a) we have a correlation between an operator corresponding to the 
first four eigenvalues and an operator corresponding to the second four eigenvalues 
in clusters A and B, respectively. In contrast, configuration (b) shows a correlation 
between operators corresponding to the largest eigenvalue only and configuration (c) 
a correlation between operators corresponding to the second eigenvalue only. 



combination of a; + 1 with x' requires the smallest number of pairs to be present in 
between the two clusters. The alternative combination is x with x' + 1, which requires 
an additional pair in between (see figure [s]). We could estimate their weights since 
the relative probability of an extra pair is the factor associated with increasing the 
separation by two. Since the correlations decay roughly as ~ 10~^ (see figure 10), we 
predict two orders of magnitude. Similarly, when x' — x is even, we get at mixture of 
the first and second four operators (see figure |8]). This explains the difference in the 
weights of the two operator sets. 

Thus, it turns out that for the IP correlations a cluster size of one rung would 
already have been large enough to reveal the dominant correlations. We will hence use 
as operator basis 

0^'± = 1^ ® -L (c^_^+^ ± ,^+1) (6.7a) 

6^'± = ^ (ct,a= ± ® Ix+i , (6.76) 

together with their hermitian conjugates. (The fact that our operator basis consists only 
of operators acting on a single rung implies that it would have been sufficient to use 
single-rung clusters. However, for the sake of consistency with the rest of our analysis, 
we retain two-rung clusters here, too.) 

The f-matrix based on these four operators (per cluster) is diagonal with equal 
entries for a given distance r. Its Fourier transform (see figure [9]) gives a result distinct 
from the Fourier transform for CD and 2P correlations. The dominant wave vectors are 
k = i/cp and k = it ± kp, where the latter is the product of an oscillation with k = ir 
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Figure 9. Fourier transform of the f-matrix obtained similarly as figure |5] for IP 
correlations based on the four operators per cluster for (a) a filling of = 0.248 and 
(b) a filling oi v = 0.286. We find peaks at about k = ±kp and k = ±kp + n (dashed 
black lines). 



and an oscillation with k = ±k-p. In total we have an oscillation in the correlations of 
the form (1 + (— l)^)e^*'^^^, i.e. an oscillation with k = ifep, and every second term 
being close to zero. The dominant wave vector /c = i/cp i s consistent with the usual 
behaviour of IP Green's functions. 

The reason for every second term being essentially zero is that the dominant hopping 
in the system, the correlated hopping, always changes the position of a particle by two 
rungs, so every second position is omitted. The small but finite value for hopping onto 
intermediate rungs is related to the finite ty/tc = 10~^ that we use. It results in a second 
oscillation at = ±kp located at intermediate rungs, whose relative strength compared 
to the dominant one is about 10~^, which is consistent with the ratio ty/tc that we used 



(see figure 10). 



We fit the one independent f-matrix element f^'^ to an exponential decay of the form 



y^g-r-Z-ri ^ggg figure 10), but apart from this we were not able to fit the exact functional 
dependence on r, especially the oscillations with k = ±kp. The reason for this is the 
existence of two oscillations where one is zero on every second rung, and that the data 
range for which reasonable IP correlations are still present is too small and thus makes 
it susceptible to numerical noise. This can be seen already in the Fourier spectrum, 
where we find relatively broad peaks, as a result of the infiuence of the exponential 
envelope and the relatively short distance range available. 

6.3. Two-particle correlations 

The operator subspace for 2P (AA^ = 2), in a cluster including two rungs has the 
comparatively small dimension of four due to the infinite nearest-neighbour repulsion 
(see figure [2]). These are c^^xCix+i, cj^,xCt,a;+i and their hermitian conjugates. In the 
present case of dominating tc, these operators represent the creation- and annihilation- 
operators of bound pairs [2J. The operator analysis yields exactly the same four 
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Figure 10. The IP correlations for a filling of (a),(c) = 0.248 and (b),(d) 
v — 0.286. Panels (a) and (b) show the IP correlations (blue symbols) together with 
a fit of the form Ae~'^^'^^ (red line). Panels (b) and (d) show the rescaled correlator 
f^'^{r)/Ae~^^'^^ (blue symbols) for distances up to r = 20. (Larger distances are 
omitted, because for these (r) < 10^^^, which is the maximal computer precision.) 
We see a strong oscillation (green curve) and a weak oscillation (red curve). 



operators with degenerate weight for all distance regimes for both cluster A and B. The 
four operators are l/\/2 {c^^xCi,x+i ± ci^xC'\,x+i) together with their hermitian conjugates, 
and they already represent the symmetric and antisymmetric combinations of the 
operators mentioned above. 

The f-matrix (4.7) is diagonal in the basis of the four operators, with equal strength 
of correlations for a fixed distance apart from a possible sign. This may be expected, 
given the similar structure of the operators. 

As for the CD correlations (AA^ = 0), we apply a Fourier transform on the f-matrix 
(see figure 11) to identify the dominant wave vectors. Again, we find two spectra of 
similar form but shifted by vr with respect to each other. Consequently we redefine 
to e^'^^f'^, the part of the f-matrix belonging to the symmetric operators. Thus, we 
obtain one leading peak at k = and sub- leading peaks at = 2k-p- Given the similar 
structure of the Fourier spectrum to that of the CD correlations, we fit the elements of 
the f-matrix to the form (6.3), but now expect 7' < 7 from the relative sharpness of the 
peaks. Already at the level of the f-matrix elements we find an overall leading decay 
with residual oscillations, whose relative magnitude becomes smaller at large distances 
(since 7' < 7). Since all matrix elements are the same after redefining /"'', it is sufficient 
to fit l/^'^l for a given /i, which will have dominant fc- vectors k = and k = ±2kp. The 
fit has errors of less than 5% throughout, with results as shown in figure 12 The overall 
behaviour is very similar to the one already found from the r.m.s. net correlations of 



Numerical results 



27 




Figure 11. Fourier transform of the f-matrix for 2P correlations based on the operators 
chosen from the four-dimensional operator space for (a) a filling of — 0.248 and (b) 
a filling oi ly — 0.286. For a detailed description see figure [5] 




Figure 12. Fitting the 2P correlations to the form in (6.3 1 for a filling of = 0.248 
and v — 0.286. The single points (blue circles and squares) are data points from the 
f-matrix and the lines (red and green) are the result of the fitting. 



this sector (see figure |4]), up to the oscillatory part from the second term in ( 6.3[ ). We 
see that the oscillations clearly decay more strongly than the actual strength |/'^''^|, in 
accord with 7' < 7. 

In contrast to the CD correlations (see figure 6.1.2), for the 2P correlations we do 
not find correlations which oscillate with phases shifted by A0 = ±7r/2 . This may come 
from the fact that clusters with the size of two rungs have the minimal possible size to 
capture 2P correlations. The corresponding operator space has dimension four and the 
four possible operators are very similar in structure. We expect that for larger clusters 
and hence a larger operator space, we would find correlations which also oscillate out of 
phase such that their oscillations cancel in the r.m.s. net correlations , in accord with 
(pl). 
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Table 3. Comparison of the power-law exponents, which we extracted from our 
numerical data, with those predicted in 





CD 




2P 


filling 




70 


W. 72 


0.248 


1.13 


1.45 


0.5 0.95 


0.286 


1.04 


1.33 


0.5 1.11 



7. Comparison to previous results 

We are now ready to compare our CDM-based results with those obtained in [2] by 
Cheong and Henley (CH) from fitting simple correlation functions. The latter were 
computed exactly in [2] for accessible separations after mapping the large model onto 
a hard-core bosonic system, but the functional forms of the r dependencies were inferred 
from a purely numerical fitting procedure. 



Overall, our results for the Hamiltonian (2.1) in the strongly correlated hopping 
regime agree with [2], in that (i) 2P correlations and CD correlations show power-law 
behaviour, (ii) the 2P correlations dominate at large distances for the fillings we were 
investigating, (iii) IP correlations are exponentially decaying and are negligible over all 
but very short distances, and (iv) the dominating fc-vector, for either 2P or CD sectors, 
is 2k-F. 



However, the power-law exponents obtained from fitting f- matrix elements to (4.10) 
and summarized in table |3| clearly deviate from the results in [2] by CH. For the CD 
correlations, in [2] the dependence of 70 on the filling u was given by the exponent 
7^^ = I + 2(2 " from which our results deviate (see figure |4] a,b) by about 25%. 
Nevertheless, our results for 70 agree quahtatively with this prediction, in that we also 
find 7o to decrease linearly with increasing filling. 

The 2P correlations deviate more strongly. For the dominant 2P correlations, CH 
predicted a constant power-law exponent of 7^^ = ^ independent of filling, coming 
from a universal correlation exponent for a chain of tightly-bound spinless fermion pairs 
[8]. In contrast, we obtain a larger exponent (see figure |4] a,b) for given fillings. Our 
result for 72 linearly decreases as the filling gets smaller and appears to approach | 
only in the limit u 0. We also explicitly calculated the same correlation function 
as investigated in |2] but found a stronger decay than the suggested there. We do 
not know whether the deviation is an artifact of the boundaries of our finite system, 
or whether the mapping used in [2] to a set of hardcore bosons might have omitted an 
important contribution. 

Moreover, it may be noted that by extrapolating the exponents in a linear fashion 
towards large fillings {u —>■ |), it appears that for fillings larger than ~ 0.35 eventually 



the CD correlations dominate over 2P correlations (see figure 13). This conclusion 
has also been found in [9] which similarly addresses diatomic real space pairing in the 
context of superconductivity. Their discussion, however, is not specifically constrained 
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Figure 13. The power-law exponents for CD correlations (70, blue symbols) and 2P 
correlations (72, red symbols) obtained from the r.m.s. net correlations for several 
fillings V. We used chain lengths of A^ = 100 (circles), A^ = 150 (crosses), and 
N = 200 (triangles). The dashed blue and red lines are linear fits to our numerical 
data for 70 and 72, respectively. The solid blue and red lines show the corresponding 
predictions of Cheong and Henley P]. For the 2P correlations, our data implies a 
linear j/-dependence going from ^ for = to | for i' = ^. This crossover from ^ to 
I is predicted by Cheong and Henley as a sub-leading contribution, without giving an 
explicit functional dependence on v. The two linear i/-dependencies imply that for large 
fillings CD correlations should become dominant over 2P correlations. Unfortunately, 
we do not have been able to obtain reliable data in that regime, because the r.m.s. 
net correlations showed strong oscillations here, contrary to our expectations from 
section 11.21 



to one- dimensional systems, and one may wonder how the specific choice of parameters 
compare. 

As the fiUing approaches 0.5 in an excluded-fermion chain, it is appropriate to think 
about the degrees of freedom as impurity states or holes in the crystalline matrix of pairs 
[9] . Then the natural length scale is the spacing between holes. The longer that spacing 
gets (it diverges as u ^ 0.5), the larger also the system under investigation must be in 
order to reach the asymptotic limit. In other words, to see proper scaling behavior in a 
uniform way, the system size should increase proportional to 1/(0.5 — i^)- In our case the 



data became unreliable for u > 0.4 (see figure 13). On the other hand, for certain fillings 



z/ < 0.4, we calculated the power-law exponents for CD and 2P correlations for ladders 



of length = 150 and = 200 (this data is also included in figure 13) and did not 



find different behaviour compared to out original data for ladders of length A^ = 100. 
8. Conclusions 



Summarizing, we found that the CDM is a useful tool to detect dominant correlations 
in a quantum lattice system. Starting from a ground state calculated with DMRG, we 
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extracted all the important correlations present in our model system. We developed 
a method which, first, determines the distance-independent operators on each cluster 
that carry the dominant correlations of the system, and second, encodes the distance- 
dependence of the correlations in the f-matrix. The latter is then analyzed in terms of 
decaying and oscillatory terms to extract the long-range behaviour of the correlations. 

We saw that the size of the clusters A and S is a hmitation of the method as it 
constrains the analysis to local operators. For some kind of correlations, however, larger 
clusters are needed to capture the relevant physics. This is not too easily implemented 
as it requires significantly more resources. As a possible alternative and as an outlook 
for possible future work, one may think of using a different cluster structure: one cluster 
as before and one "super-cluster" representing a larger continuous part of the system 
including one boundary. As MPS introduces, for each site, effective left and right 
Hilbert spaces describing the part of the chain to the left and to the right of that site, 
the description of such a super-cluster should be straightforward. The resulting effective 
density matrix describing a large part of the system can be calculated accordingly. 

Overall, DMRG is a suitable method to calculate the CDM. The latter is easily 
and efficiently calculated within the framework of the MPS. The explicit breaking of (i) 
translational invariance by using finite system DMRG and (ii) a discrete symmetry of 
the model, lead us to develop certain strategies to restore these broken symmetries. The 
smoothing of the boundaries can still be further optimized, or be replaced by periodic 
boundary conditions. However, we do not expect that this will have significant influence 
on the conclusions drawn. 
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A. The veiriational matrix product state approach 

This appendix offers a tutorial introduction to the variational formulation of DMRG for 
finding the ground state of a one-dimensional quantum lattice model, , based on matrix 
product states (MPS). It also explains how this approach can be used to efficiently 
calculate the GDM. We point out all the important properties of the MPS and explain 
how to perform basic quantum calculations such as evaluating scalar products and 
expectation values, as well as determining the action of local operators on the MPS and 
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constructing a reduced density matrix. We explain how a given MPS can be optimized 
in an iterative fashion to find an excellent approximation for the global ground state. 
We also indicate briefiy how the efficiency of the method can be enhanced by using 
Abelian symmetries. 

We would like to emphasize that we make no attempt below at a historical overview 
of the DMRG approach, or at a complete set of references, since numerous detailed 
expositions of this approach already exist in the literature (see the excellent review by 
U. Schollwock [1]). Our aim is much more modest, namely to describe the strategy 
implemented in our code in enough detail to be understandable for interested non- 
experts. 

A.l. Introduction 

Quantum many-body systems deal with very large Hilbert spaces even for relatively 
small system sizes. For example, a one- dimensional quantum chain of spin ^ 
particles forms a Hilbert space of dimension 2^, which is exponential in system size. For 
quantum lattice models in ID a very efficient numerical method is the density matrix 
renormalization group (DMRG), introduced by Steven R. White [3]. The problem of 
large Hilbert space dimension is avoided by an efficient description of the ground state, 
which discards those parts of the Hilbert space which have negligible weight in the 
ground state. In this manner the state space dimension of the effective description 
becomes tractable, and it has been shown that this produces excellent results in many 
quasi one-dimensional systems. 

The algebraic structure of the ground state for one-dimensional systems calculated 
with DMRG is described in terms of matrix product states (MPS) [IDl HB [HI El 113] ■ The 
origin of this MPS structure can be understood as follows (a detailed description will 
follow later): pick any specific site of the quantum lattice model, say site k, representing 
a local degree of freedom whose possible values are labeled by an index ak (e.g., for a 
chain of spinless fermions, 0"^ = or 1 would represent an empty or occupied site). Any 
many-body state lip) of the full chain can be expressed in the form 

1^) = E ^trl h) Wk) \n) , (A.l) 

where \lk) and \rk) are sets of states (say Ni and Nr in number) describing the parts of 
the chain to the left and right of current site k, respectively, and for each a^, A^'^''^ is a 
matrix with matrix elements dimension A^; x Nr. Since such a description is 

possible for any site k, the state {ip) can be specified in terms of the set of all matrices 
y[[<^k]^ resulting in a matrix product state of the form 

\^)= J2 {A^''<..A^''^^),^rJ^^)---M ■ (A-2) 

0-l...(T]V 

One may now seek to minimize the ground state energy within the space of all MPS, 
treating the matrix elements of the A-matrices as variational parameters to minimize the 
expectation value {ipl H lip). If this is done by sequentially stepping through all matrices 
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in the MPS and optimizing one matrix at a time (while keeping the other matrices 
fixed), the resulting procedure is equivalent to a strictly variational minimization of the 



ground state energy within the space of all MPS of the form (A. 2) [5l [TDl ttH [121 |T3]. If 
instead the optimization is performed for two adjacent matrices at a time, the resulting 
(quasi-variational) procedure is equivalent to White's original formulation of DMRG 
[ini im HSl [I3]. The MPS based formulation of this strategy has proven to be very 
enlightening and fruitful, in particular also in conjunction with concepts from quantum 
information theory [5]. 

In general, such an approach works for both bosonic and fermionic systems. 
However, to be efficient the method needs a local Hilbert space with finite and small 
dimension, limiting its applicability to cases where the local Hilbert space is finite 
dimensional a priori (e.g. fermions or hard-core bosons) or effectively reduced to a finite 
dimension, e.g. by interactions. For example, such a reduction is possible if there is a 
large repulsion between bosons on the same site such that only a few states with small 
occupation number will actually take part in the ground state. For fermions, on the 
other hand, the fermionic sign must be properly taken care of. The anti-commutation 
rules of fermionic creation and annihilation operators causes the action of an operator 
on a single site to be non-local because the occupations of the other sites have to be 
accounted for. To simplify the problem, a Jordan- Wigner transformation [13] can be 
used to transform fermionic creation and annihilation operators to new operators that 
obey bosonic commutation relations for any two operators referring to different sites. 
This greatly simplifies the numerical treatment of these operators as fermionic signs can 
be (almost) ignored. 

Before outlining in more detail the above-mentioned optimization scheme for 
determining the ground state (see section A. 3), we present in section A. 2 various 



technical ingredients needed when working with MPS. 
A . 2. Matrix product states 

A. 2.1. Construction of matrix product states We consider a chain with open boundary 
conditions consisting of equal sites with a local Hilbert space dimension of d. A state 
lip) is described by 

1^) = ^'^i'-,<^iv ki) • • • Wn) , (A.3) 

(Tl...(Tff 

where = 1,. . . ,d labels the local basis states of site i. In general, the size of the 
coefficient space ip scales with 0{d^). This can be rewritten in a matrix decomposition 



of the form (A. 2) with a set of times d matrices A''^*' (see section A. 2. 3 for details). 
Formally, this decomposition has two open indices, namely the first index of A^'^^J and 
the second index of A'""^', as A^'^i] and A^'^^l are not multiplied onto a matrix to the left 
and to the right, respectively. For periodic boundary conditions these two indices would 
be connected by a trace over the matrix decomposition, giving a scalar. In the case of 
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open boundary conditions, the two indices range only over one value (see section A.2.3), 
i.e. the matrix decomposition is a 1 x 1 matrix which is a scalar. 

If these A-matrices are sufficiently large this decomposition is formally exact, but 
since that would require A-matrices of exponentially large size, such an exact description 
is of academic interest only. The reason why the A-matrices are introduced is that they 
offer a very intuitive strategy for reducing the numerical resources needed to describe a 
given quantum state. This strategy involves limiting the dimensions of these matrices by 
systematically using singular-value decomposition and retaining only the set of largest 
singular values. The A-matrices can be chosen much smaller while still giving a very 
good approximation of the state \ip)- 



Selecting a certain site k, the state can be rewritten in the form (A.l). The 
effective 'left' basis \lk) = Y.ai...ak-i ^^"^"^ ■ ■ ■ Wi) ■ ■ ■ Wk^i) describes the sites j = 

1, . . . , k — 1, the effective 'right' basis \rk) similarly describes the sites j = k + 1, . . . ,N. 
Site k is called the current site, as the description of the state makes explicit only the 



A-matrix of this site (see figure Al ) 






[k- 1) — j-\^J—]^ — + 1) 

h) IcTfc) kfc) 

Figure Al. Current site witii effective basis sets. 



So far (A. 3) and (A.l ) are equivalent, but now we have a representation of the state 
which allows a convenient truncation of the total Hilbert space, used for the description 
of a MPS. For example, if we introduce a parameter D and truncate all effective Hilbert 
spaces of all sites to the dimension D, each A^'^'^l-matrix has at most the dimension DxD. 
This reduces the resources used to describe a state from 0{d^) for the full many-body 
Hilbert space down to 0{N D'^d). This is linear in the system size, assuming that 
the size required for D to accurately describe the state grows significantly slower than 
linearly in A^. This, in fact, turns out to be the case for ground state calculations [15]. 
Details of this truncation procedure and estimates of the resulting error are described 
in section IA.2.5[ 

A. 2. 2. Global view and local view Matrix product states can be viewed in two 
alternative ways: a global view and a local view. Both views are equivalent and both 



have their applications. In the global view the state is expressed as in (A. 2), i.e. the 
effective Hilbert spaces have been used 'only' to reduce resources. The state is stored 
in the A-matrices, but the effective basis sets will be contracted out. This perception 
has to be handled very careful, because contracting out the effective basis sets leads to 



higher costs in resources! In the local view the state is expressed as in (A.l). It is called 
local because there is one special site, the current site, and all other sites are combined 
in effective orthonormalized basis sets. Usually, the local view is used iteratively for 
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every site. In this perception, we need effective descriptions of operators contributing 



to the Hamiltonian acting on other sites than the current site (see section A. 2. 8). 



A.2.3. Details of the A-matrices The A-matrices have some useful properties that 
hold independently of the truncation scheme used to limit the effective Hilbert spaces. 
First of all, we notice that by construction dim{T-C'^-^) = dim{T-C'''=), otherwise the matrix 



products in (A. 2) would be ill defined. Based on this, we can find another interpretation 
of the y4-matrices in the local view. The part of the chain to the left of site k (where k 
is far from the ends for simplicity) is described by the effective basis \lk), which is built 
of truncated A-matrices: 



o"lviCrfc_i 



EE E ('""■'■■■^""')„_.K>---K-)ifc;:ik-.> 



S 



|/fc-l> 

= E ^tik\^^-^)M ■ (A-4) 

The y4[°"''-i]-matrix maps the effective left basis together with the local \<Jk-i) basis 
onto the effective left basis The same argument applied on the effective right basis 
of site k leads to the transformation of \rk+i) and |crfe_|_i) onto \rk) via the y4[°"'=+il-matrix: 

So far, this may be any transformation, but in order to deal with properly orthonormal 
basis sets, we may impose unitarity on the transformation (see below). 

The A-matrices towards the ends of the chain have to be discussed separately. The 
use of open boundary conditions implies that we have a 1-dimensional effective state 
space to the left of site one and the right of site A^, respectively, both representing the 
empty state. This implies that dim(7i'^) = 1 = dim('W^). Moving inwards from the 
ends of the chain, the effective Hilbert spaces acquire dimension d^jd"^,... until they 
become larger than D and need to be truncated. Correspondingly, the dimension of 
matrix At""*] is D^-i x D^, where Di~ = min((i^, d^~'', D). There is no truncation needed 
if dim(7^''-) *d = dimiW") or dim(H'''=) *d = dim(7^'*). In these cases we simply choose 
Aikcrk)rk = 1 and Ai^^r.a,) = 1, respcctivcly. 



Summarizing, the A-matrices have two functions. If site i is the current site in ( A.l ), 
the A ["^'l -matrices represent the state, i.e. its coefficients specify the linear combination 
of basis states \ak) and jr^). On the other hand, if not the current site, the A- 
matrices are used as a mapping to build the effective orthonormal basis for the current 
site, as we describe next: 



Orthonormal basis sets In the local view, the whole system is described by the A- 
matrices of the current site k in the effective left basis, the effective right basis, and the 
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local basis of site k. A priori, the basis states form an orthonormal set only for the local 
basis set, but we may ask for the effective basis sets |/) and |r) [|] to be orthonormal, 
too, i.e. require them to obey: 

(r'|r) = 6r'r ■ (A.6) 



This immediately implies the following condition on the ^''^^l-matrices, using (A. 4) and 



(A. 5) (for a derivation, see section A. 5.1): 

=1 for J > A;. 



(A.7) 



The orthonormality (A.6) for both the left- and right basis states holds only for the 



current site. For the other sites there is always only one orthonormal effective basis. 



Graphical representation Matrix product states can be depicted in a convenient 



graphical representation (see figure A2 ). In this representation, A-matrices are displayed 
as boxes and A^'^''^ is replaced by Af^ for brevity. Indices correspond to links from the 
boxes. The left link connects to the effective left basis, the right link to the right one, 
and the link at the bottom to the local basis. Sometimes indices are explicitly written on 
the links to emphasize the structure of the sketch. Connected links denote a summation 
over the indices (also called contraction) of the corresponding A t*^] -matrices. At the 
boundaries of the chain, a cross is used to indicate the vacuum state. 



(a) 



(b) 



x-|^3~~l^3 



0"2 



cn-i Cat 



'^1^* i/fc) = X— 1^3 — "11*3" I'^'t) = "E^ — 



0"i 



0-1 



{Tat 



Figure A2. Graphical representation of a matrix product state in the (a) global view 
and (b) local view. 



A. 2. 4- Orthonormalization of effective basis states We now describe how an arbitrary 
MPS state can be rewritten into a form where its local view with respect to a given site 
has orthonormal left- and right basis states. It should be emphasized that this really just 
amounts to a reshuffling of information among the state's A-matrices without changing 
the state itself, by exploiting the freedom that we always can insert any X~^X = 1 at 
any position in the matrix product state without altering it. 

I From now on the index k is only displayed when several sites are involved. For the current site or in 
the case when only one A-matrix is considered the index will be dropped. 
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Assume site k to be the current site and assume that it has an orthonormal left basis 
(the latter is automatically fulfilled for k = 1). We need a procedure to ensure that, when 
the current site is switched to site k+1, this site, too, will have an orthonormal left basis. 



(This is required for the orthonormality properties used in the proof in section A. 5.1 A 
similar procedure can be used to ensure that site k — 1 has an orthonormal right basis 
provided k has su ch a b asis.) For this purpose we use the singular value decomposition 
(SVD, see section A. 5. 2) for which we have to rewrite fusing the indices Ik and 

(7k- 



A 



m,n 



E 



mrk 



(A.8) 



where m, n and have the same index range (see figure A3). Specifically, u fulfills 

1=U^U= u\i^ak),-m'U{lkak),m , (A.9) 



{IkO-k) 



which is equivalent to the orthonormality condition (A. 7) for the A^'^'^^-matrices. 
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Figure A3. Singular value decomposition of the A-matrices 
As u replaces A^"'''^ and sv"^ is contracted onto A^"'^+'^\ this leaves the overall state 



unchanged (for a graphical depiction see figure A4): 

^ ^ - ^IkTk^lk+irk+l ~ Z^^'lkmV^'J Jmrk Ik+lVk+l 

{rk=lk+i) (rk=h+l) ^ 

= (sv^Ak+i) = . (A. 10) 
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Figure A4. Rearrangement of the A-matrices to switch the current site from site k 
to fc+ 1. 

Site k + 1 now has an orthonormal effective left basis. A similar procedure works for 



the effective right basis, see figure |A5[ To obtain an orthonormal effective left basis for 

4-1 ■ 
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Figure A5. Orthonormal effective right basis for site k — 1. 
the current site k, we start with the first site, update A^'^^^ and A^'^^\ move to the next 
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site, update A^°'^^ and A^'^^\ and so on until site k — 1. For an orthonormal effective right 
basis, we start from site and apply an analogous procedure in the other direction. 

If the state is in the local description of site k with orthonormal basis sets 
\h), Wk) and \rk), it is now very easy to change the current site to site ± 1, with 
corresponding new orthonormal basis sets \lk±i), \o'k±i), \rk±i)- Suppose we want to 
change the current site from site k to site fc + l. Following the procedure described 
above, site fc + l already has an orthonormal right basis and all sites left of site fc fulfill 
the orthonormality condition. All that is left to do, is to update site fc and fc + 1 to 
obtain an orthonormal left basis for site fc + l. This is called a switch of the current site 
from site fc to fc + 1. The switch from site fc to site fc — 1 is done analogously. 

A. 2. 5. Hilbert space truncation A central ingredient in the variational optimization of 
the ground state (see section A.3.1| below) is the truncation of the effective Hilbert spaces 



associated with a given A-matrix. The strategy for truncating the effective Hilbert 
spaces is completely analogous to the original DMRG formulation [TT]. The DMRG 
truncation scheme is based on discarding that part of the Hilbert space on which a 
certain density matrix has sufficiently small weight. There are two ways to obtain an 
appropriate reduced density matrix: two-site DMRG [31 H] and one-site DMRG The 
crucial difference between the two is that one-site DMRG is strictly variational in the 
sense that the energy is monotonically decreasing with each step,, whereas in two-site 
DMRG the energy may (slightly) increase in some steps, but with the advantage that 
the cutoff dimension can be chosen dynamically in each step. 

Two-site DMRG Two-site DMRG arises when variationally optimizing two sites at a 
time. We consider two current sites, say fc and fc + l, and we may choose the cutoff 
dimension site-dependent: D Dk = dimilH}*') . Following section A. 2. 4, we assume 



site fc to have an orthonormal left basis and site fc + 1 to have an orthonormal right 



basis. After contracting the indices connecting A^"'''^ and A^'^^+^^ (see figure A6), the 
state is described by ^/^r^+T^^' "'■^ ^^^^ description we may optimize the ground state 
locally by variationally minimizing the ground state energy with respect to ^1^*^^^^'' 
(see section A. 3.1). Afterwards, we need to decompose into ^''^'=1 and y4[°"'=+il 



again. This can be accomplished via singular value decomposition (see section A. 5. 2) by 



fusing the indices /fc,o"fc — >■ {lk<^k) and rk+i,<^k+i ('^fc+iO'fc+i) (see figure A6) to obtain 

I kfe+i] 



AtT+r' = where i = 1 . . .mm{dDk,dDk+2). Using the column 



unitarity of u and the row unitarity of v' (see section A. 5. 2), we rewrite the state as 




(^o!::;;' j h) Wk) Wk^i) ir.+i) 
J2 (^^)!r,t!' i^^+i^ I'^^+i^ 



^•fc+lO'fe+l 

^ — 

\ri) 
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= 1^5.101^.) , (A.ll) 

i 

where the new set of basis states and |fj) is orthonormal with = and 

i^i'l'f^i) = Si'i- This representation of the state may be seen as residing on the bond 
between k and k + 1, with effective orthonormal basis sets for the parts of the system to 
the left and right of the bond. Reduced density matrices for these parts of the system, 
obtained by tracing out the respective complementary part, have the form: 

i i 

The standard DMRG truncation scheme amounts to truncating p'-^l and pt^' according 
to their singular values Sj. We could either keep all singular values greater than a 
certain cutoff, thereby specifying a value for Dk+i between 1 and min ((iZ)fc, or 
alternatively choose Dj. = D iohe site-independent for simplicity. This step makes the 
method not strictly variational, since we discard some part of the Hilbert space which 
could increase the energy. It turns out that this potential increase of energy is negligible 
in practice. We can obtain a measure for the information lost due to truncation by using 
the von Neumann entropy S" = — tr (plnp), given by 

i>D 

where ^ = 1 due to the normalization of 
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Figure A6. Procedure for site update within two-site DMRG. The grey hue under 
the s indicates that s is the diagonal matrix of singular values. 



One-site DMRG One-site DMRG arises when variationally optimizing one site at a 
time. In contrast to two-site DMRG, one-site DMRG does not easily allow for dynamical 
truncation during the calculation. (It is possible in principle to implement the latter, 
but if one decides to use dynamical truncation, it would be advisable to do so using two- 
site DMRG.) The truncation is fixed by the initial choice of but it is still possible to 
determine an estimate on the error of this truncation by analyzing the reduced density 
matrix. Starting from an expression for the full density matrix in the local view (current 
site k with orthonormal effective basis sets) 

p = 1^) (^1 = ( E^l? 10 1^) 1^) ) ( E A^v* (^'1 (^'1 ) 

V Ira J Xl'r'a' / 

= E ^I?^'* 10 (n k) \r) {r'\ , (A.14) 

Iral'r'a' 
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we trace out the effective right basis and obtain a reduced density matrix for the current 
site and the left part of the system: 



fc+ij 



E 10 (n w) w\ ■ (A.15) 



p 

Iral'cr' 

This reduced density matrix carries the label Ik+i because it corresponds precisely to 
the density matrix |/a;+i) (^fc+i |- So if we switch the current site from site k to site k + 1, 
we can check the error of the truncation of Ti^'-'+i . Fusing the indices I and a, we obtain 

Iral'a' Iral'a' 
lal'a' 

We do not need to diagonalize the coefficient matrix AA'' to obtain the largest weights 
in the density matrix, because we get its eigenvalues as a byproduct of the following 
manipulations anyway [4J. To switch the current site we need to apply a singular 



value decomposition (see section A. 2. 4) and obtain A = usv'^ (this is not the usual A- 
matrix, but the index- fused form). This directly yields AA^ = usv'^vsv) = us'^u\ which 
corresponds to the diagonalization of implying that the weights of the density 

matrix are equal to s^. Of course this works also for the right effective basis. With such 
an expression, we can check whether the effective Hilbert space dimension D of is 
too small or not. For example, we could ask for the smallest singular value sd to be at 
least n orders of magnitude smaller than the largest one Si, i.e. the respective weights 
in the density matrix would be 2n orders of magnitude apart. If the singular values do 
not decrease that rapidly, we have to choose a greater D. 

A. 2. 6. Scalar product The scalar product of two states and \ip') is one of the 
simplest operations we can perform with matrix product states. It is calculated 
most conveniently in the global view because then we do not need to care about 
orthonormalization of the A-matrices: 

a'^...a'j^ CT1...CTJV 

= E (^''"'' • • • ^''""^T (^^"'' • • • ^^"^0 ' (A- 17) 

(Ti...crjV 

using the ort honor mality of the local basis (cr[,|o"i) = ^fc/^o-^o-fc- In principle the order 
in which these contractions are carried out is irrelevant, but in practice it is possible 
to choose an order in which this summation over the full Hilbert space is carried out 
very efficiently by exploiting the one-dimensional structure of the matrix product state 



(see figure A7 for a graphical explanation). For details on the numerical costs, see 
section A. 5. 3, In method (a), after contracting all A-matrices of and we have 
to perform a contraction over the full Hilbert space, i.e. a 1 x c?^ matrix is multiplied 
with a c?^ X 1 matrix. This contraction is of order O [d^) , which is completely unfeasible 
for practical purposes. In method (b) the most 'expensive' contraction is in the middle 
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of the chain, say at site k, and it is of order O (dD^). Here the A-matrices are viewed 
as three-index objects Ai^^^^^ with dimension D x D x d. All sites left of site k are 
represented hj a D x D matrix, say L^f. Contracting this with the matrix at site k 

yields the object L'Mij.rfc<Tj. , which has dimensions D x D x d, and since the sum 
contains D terms, the overall cost is 0{dD^). Thus, in practice, method (b) is rather 
efficient and renders such calculations feasible in practice. 
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Figure A7. Scalar product, computed in two different orders, (a) First all ^-matrices 
of and are contracted and then contraction over the local indices is carried 
out. b) First, for site one, we contract over the local indices of Ai and A'^. Then we 
contract over the effective index between Ai and A2 and afterwards over the indices 
between the resulting object and {A'2)* ■ Proceeding over the whole chain yields the 
scalar product. 



Partial product Sometimes it is required to calculate a product over only a part of the 
matrix product state. This is done the same way as the scalar product 

(p[^'=l)^ ^, = (^["^1 . . . ^["'=-^1)* (A^'i] . . . ^["'=-^1)^ , (A.18) 



(7l...O-fc_i 



(A. 19) 



o-fc + i...crjv 



p[kk'] 



V . . . ^["^'-il)*, ,, (^["''■+11 . . . At'^'-il) , . (A.20) 

''fe 'fe' "^k '■y 



o"fc + l..-o-fc'_l 



Notice that P^^"^ and P^^'=^ are matrices in the indices Ik and r^, respectively (see 



figure A8). In fact, they correspond to the overlap matrices (I'^lh) and {r'k\rk), 
respectively. 
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Figure A8. Partial products associated with site k. 



A. 2.1. Reduced density matrix The pure density matrix given by the matrix product 
state \ip) is defined as p = {ip]. To describe only a part of the system, we need to 
calculate the reduced density matrix. Let / be a set of sites and as = {^kei} a fused 
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index for their local states. Tracing out all other sites with combined index = {o"a,.^/} 
we obtain 

pj= Yl K< (^'''^' • • • ^^'^"0* (^^'"'^ • • • ^'''"'^ '^^^ • ^^-^^^ 

a^...aj^cT[...a'j^ 

This is a completely general expression, but in the cases where / = {k} or I = {k,k'} 



it reduces to (see figure A9 ) 

P{k} = P^"-'^ (a'"'^' ® A'"^^*) P^""'^ Wk) (^fcl , (A.22) 

p^,,,j = p[^^] (^Al-^^l ® AKI*) pl'^'^'l (^A^-^'] ® aK']*) pl«.'] |afe) {cr',\ {a',,\ . (A.23) 

A similar strategy can be used to calculate the density matrices needed for the main 
text, by contracting out the a^s for all sites except those involved in the clusters A, B 
or y4 U -B. In fact, (A.23) gives p^^^ for two clusters of size one at sites k and k' . 
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Figure A9. Reduced density matrix (a) for site k and (b) P{kk>} for sites k and 
k\ where k < I < k' . 



A. 2. 8. Operators in an effective basis Let k be the current site with orthonormal 
effective basis sets \lk) and \rk). Consider an operator B, which acts on the local basis 
of site k — 1 only, with matrix elements B^-'^^^^^^^ = B \ak-i). We call this the 

{k — 1) -local-representation of B. To represent B in the effective left basis of site k, 
called the k -left-representation of 5, we use the transformation properties of A^"''-^^ 



(see figure AlO) 



I 'fe-iffc-i 



'fc-i"fc-i 



(A.24) 

'fc-lf fe_l''"fc-l 

where the only condition to derive these results, was that site k — 1 has an orthonormal 
effective left basis. Similarly, if the (fc — l)-left-representation of an operator C is known, 

(A.25) 



its /c-left-representation can be obtained via (see figure AlO) 



Equation (A.24) and (A.25) can be used iteratively to transcribe the z-local- 



representation of B into its A;- left-representation for any k > i (see figure All). This 



Variational matrix product state approach 



42 



(a) 



A, 



B 



B 



(b) 



c 



h-i 



c 



Figure AlO. The fc-left-representation of (a) the operator _B, obtained from its 
(k — l)-local-representation and (b) the operator C, obtained from its (k — l)-left- 
representation. 



Be 



Figure All. Iterative calculation of the /c-left-description of an operator _B, given in 



the i-local-description, by (A. 24 1 and (A. 25 1 for any k> i 



reasoning also applies to the right site of site k and so it is possible to obtain a description 
of any local operator on any site. 

To obtain a description of a pair of local operators acting on different sites, we 
have to transcribe them step by step. Let site k be the current site with orthonormal 
effective basis sets and -B, C two operators acting locally on site i and j respectively 
{i < j < k). First we obtain the j-left-representation of B, namely Bi'i-, as described 
above. Then both operators are transformed together into the (j + l)-left-representation 



(see figure A12) 



E 



(A.26) 



which in turn can be transformed iteratively into the desired fc-left-representation of the 
operators B and C. 







A- 


1 1 










a 


; 




B 




c 




BC 








a 


/ 

; 








^; 


— 





'3+1 



Figure A12. The (j + l)-left-representation of the operators C, given in the j-local- 
representation, and B, given in the j-left-representation. 



A. 2. 9. Local operators acting on {tp) Any combination of operators can be calculated 
directly in the global view or in the local view via the effective descriptions introduced 
in the previous section. 
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Global view The operators, known in the local basis of the site they are acting on, are 
contracted directly with the corresponding A-matrix. For example, the formula for a 



nearest neighbour hopping term clck+i (see figure A13) reads as 



. . . \ak-i) \ak+2) ■ ■ ■ Wn) ■ (A.27) 



Local view Let k be the current site with orthonormal effective basis sets. If we want 
to evaluate operators acting on other sites than the current site k, we need an effective 
description of these operators in one of the effective basis sets of site k to contract these 
operators with the A-matrix of the current site. For example, to calculate the action of 
the nearest neighbour hopping term clck+i on = Aj^*^ \l) \ak) \r), we need (c|,)o-' 



and {ck+i)r'r to obtain (see figure A13) 



E ( E (4)^,^^ ) ( E ) A^^ 10 K) \r') ■ (A.28) 



(a) x-f^; 



Ao A 



Ak+i\ — piV-1 

r zi 



Cfc+l 



a2^ (b) -^4 



Ck + l 



0\ 



0-2 



Figure A13. The nearest neighbour hopping term c^Cfc+i acting on j?/;) in (a) the 
global view and (b) the local view. 



A.2.10. Expectation values Expectation values are merely the scalar product between 
the state with itself including the action of an operator and can be easily worked out in 



both the global and the local view (see figure A14). Since both methods are equivalent, 
the local variant is much more efficient as it involves much less matrix multiplications. 
However, it requires careful orthonormalization of the remainder of the A-matrices. The 
iterative scheme, introduced in section A. 3 allows for that and works in the local picture. 
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Figure A14. The expectation value of the nearest neighbour hopping c^c/j^i in (a) 
the global view and (b) the local view. 
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A. 3. Variational optimization scheme 



The basic techniques introduced in the previous sections are the building blocks for 
DMRG sweeps, an iterative scheme to determine the ground state in the usual DMRG 
sense. This scheme starts at some site as current site, for example the first site where 
truncation occurs, and minimizes the energy of {ip) with respect to that site. Afterwards 
the current site is shifted to the next site, and the energy of lip) with respect to that site 
is minimized. This is repeated until the last site where truncation occurs is reached and 
the direction of the switches is reversed. When the starting site is reached again, one 
sweep has been finished (see figure A15). These sweeps are repeated until \ip) converges. 



A, 



truncated ^—matrices 



. A 



Figure A15. One complete sweep. 



A. 3.1. Energy minimization of the current site In order to find the ground state of the 
system we have to minimize the energy E = {ip\H of the matrix product state {ip) 
with the constraint that the norm of \ip) must not change. Introducing A as Lagrange 
multiplier to ensure proper normalization, we arrive at the problem of determining 



Tam{{ilj\H\ 



-X{^m. (A.29) 
In the sweeping procedure introduced above, the current site is changed from one site 



to the next and the energy is minimized in each local description. Thus, we need (A.29) 



in terms of the parameters of the current site. Let us describe how to do this for the 
case of one-site DMRG, where the A-matrices are optimized one site at a time. (The 
procedure for two-site DMRG is entirely analogous, except that it involves combining 
A-matrices of two neighboring sites by fusing their indices to obtain a combined two-site 



y4-matrix, see section A. 2. 5 ) Inserting (A.l) into (A.29) yields (see figure A16) 



\lral' r' a' Ira / 



(A.30) 



where Hiir'a'ira = {^'\ {(^'l {^'\ H \l) \a) \r) is the Hamiltonian expressed in the two 
orthonormal effective basis sets and the local basis of the current site. 

The multidimensional minimization problem ( |A.29 ) has been transformed to a local 
minimization problem where one ^-matrix (or two) is optimized at a time and all others 
are kept constant. Such a procedure could, in principle, cause the system to get stuck 
in a local minimum in energy, but experience shows that the procedure works well jl], 
especially in the presence of a gap. 
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Figure A16. The minimization problem expressed in the current site. 



To obtain a solution for (A. 30), we differentiate tlie equation with respect to A^^J 



(this is possible because the Hilbert space has an hermitian scalar product) and obtain 
= J2 Hvr'.'lraAf} - XA^) . (A.31) 



/'r'cr' 



The matrix elements Hi/r'a'ira may be calculated easily using the techniques introduced 
in section 



A.2 



, see section 



A. 3. 2 for details). Changing to matrix notation and replacing 



A with Eq in anticipation of its interpretation as an energy, we obtain an eigenvalue 
equation: 

HAt^ \l) \a) \r) = E,At^ \l) \a) \r) . (A.32) 

The minimization problem reduces to a local eigenvalue problem, which can be solved 
by standard techniques. The full Hilbert space of the current site has dimension dD^ 
and may become large, but it is not necessary to determine the full spectrum of H, 
since we are interested only in the ground state. The Lanczos algorithm is an effective 
algorithm to achieve exactly that. The advantage of this algorithm is that we only have 
to compute H lip), which saves much effort. The Lanczos algorithm produces as output 
the ground state eigenvalue and eigenvector. The latter gives the desired optimized 
version of the matrix Af^, which then has to be rewritten (with or without Hilbert space 
truncation, as needed) into a form that satisfies the orthonormality requirements of the 



left and right basis sets, as described in section |A.2.4 



A. 3. 2. Sweeping details Before the actual sweeping may be started we have to set 
up an initial state, prepare a current site with orthonormal effective basis sets and 
calculate effective descriptions of operators which are part of the Hamiltonian. After 
this initialization we may determine the ground state with respect to this current site 
and shift the current site to the next site. That current site again has orthonormal 



effective basis sets due to the switching procedure introduced in section |A.2.4[ but we 



also need effective representations of the operators acting in the Hamiltonian. At this 
step the structure of the matrix product state saves much effort, as most of the needed 
representations are already calculated. 



Structure of the Hamiltonian terms The Hamiltonian Hiir'a'ira, acting in the space 
spanned by the states |/), \a), |r), breaks up into several terms: 



Hl'r'a'lra 



h'l ® (^.)aV ® Ir'r + (^^l)p« ® K'a ® Ir'r + 1 



I'l 



K'a 
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+ ® ^r'r + h'l ® [H.r),,,,,, + {HL.R)nr'r,,^ , (A.33) 

where the indices denote on which parts of the system the respective term acts on (L 
and R indicate left and right of the current site, respectively, • indicates action on the 
current site). Of course, the six terms of (A.33) depend on the current site k: h'^^\ 
Hf^\ H^j'^J, hI^^ and -ff^ti?- The terms {Hl)i'i and {Hjij^',- contain all terms which 
involve only sites k' < k and k' > k, respectively. The iterative structure of the method 
directly yields the following equalities: 



H 



(fc-i) 

R 



rik) 



(A.34) 
(A.35) 



where the terms on the right hand side are meant to be expressed in the effective basis 



of the operator on the left hand site (see figure A17). 
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Figure A17. Iterative calculation of the operator H)^ 



{k+i) 



where H^l acts only on sites k' < k and Hf'^ 



that H"^^' has the form H^ i ® K 
only on site k. The calculation of Hj^ works analogously. 



The sum over i indicates 
(fc) 



Initialization First of all we need an initial matrix product state, which is most 
conveniently chosen to consist of identity transformations at the ends of the chain 
(see section A. 2. 3) and random A-matrices for the rest of the chain. We take the 
first site where Hilbert space truncation is applied as current site k and obtain an 
orthonormal effective right basis (the effective left basis is already orthonormal) using 
the orthonormalization procedure introduced in section |A.2.4 starting from site A^. 
Additionally it is convenient, while dealing with site A^, to calculate and store the 
operator H^~^^ (see (A.35)) and the effective description of all operators of site A^ 



which contribute to H. 



W 

•R 



and H^J^ in the effective right basis of site — 1 



section 



A. 2. 8). This ensures, when the sweeping procedure reaches site A^ 



(see 
that 



all necessary operators are already calculated. This is repeated from site down to 
site k + 1, and similarly for the sites k' < k in the other direction. The result of these 
initialization steps is that we have a current site k with orthonormal effective basis sets. 



effective descriptions of the Hamiltonian terms H^j^^ and H^^' and effective descriptions 



of all operators contributing to H^iJ , H^^^ and H)^Jj^. Moreover, with an appropriate 
extension to the switching procedure of section A. 2. 4, all effective descriptions for other 
current sites are available for use when needed in future sweeping steps. 



(fc) 
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Extended switching procedure The switching procedure of section A. 2.4 is apphed as 



before. Additionally, depending on the direction of the switch, ifj^^'^^^or H^^ are 



calculated and stored as well as the operators needed for the Hamiltonian (A.33). This 
extended switching ensures that for the new current site all required operators are 
calculated, if they had been for the old current site. 

Complete ground state calculation The methods introduced above make the procedure 
to determine the ground state very efficient as the global problem is mapped onto 
many local problems involving only a few terms to calculate. The iterative structure 
of the matrix product states and the effective Hamiltonian terms strongly increase the 
efficiency. A full ground state calculation consists of: 

(i) Initialization as described above 

(ii) Full sweeps from site K to site K' and back to site K, with sites K and K' the first 
and last site where the effective Hilbert spaces are truncated. 

(iii) After each sweep i the overlap {ipi-ilipi) between the state before and after the 
sweep is calculated. If the matrix product state does not change any more, stop 
the sweeping. A criterion, for example, for when to stop would be to require that 

where e is a small control parameter, typically of order 10~^°. 

Numerical costs The step with the most impact on the numerical costs of the algorithm 
is the calculation of H in the Lanczos method. This method is an iterative scheme 
using several Lanczos steps, of which usually less than 100 are needed for one ground 
state calculation. Each Lanczos step calculates H lip) exactly once. This calculation 



basically consists of elementary matrix multiplications, see section A. 5. 3 for details on 
the numerical costs of such calculations. The six terms introduced in (A.33) are not all 
equally time consuming. Most of them contain identity maps which do not need to be 
carried out and thus the term Hl,r is the most time consuming, requiring operations 
of order 0{dD'^{2D + d)). The total numerical cost for the minimization process is 

C = ATg^eep X X A^Lanczos X {dD^ {2D + d)) , (A.37) 

where N^^eep is the number of sweeps, the chain length and ALanczos the number of 
Lanczos steps. In practice the cutoff dimension is significantly higher than the local 



Hilbert space dimension d and thus (A.37) is nearly linear in d. 



A. 4- Abelian symmetries 

Matrix product states can be easily adapted to properly account for conserved quantum 
numbers, representing the global symmetries of the Hamiltonian. We will limit ourselves 
to Abelian symmetries, meaning that the irreducible representation of the symmetry 
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group is Abelian, as these are easily implemented, which is not necessarily the case for 
non-Abelian symmetries [T6] . 

An Abelian symmetry allows a quantum number Q to be attached to every state. 
The property that the symmetry is Abelian manifests itself in that this quantum number 
is strictly additive. For two states and IQ2), the quantum number of the direct 
product of these two states is given by \Qi) ® \Q2) = \Qi + Q2)- For example, if the 
Hamiltonian commutes with the number operator for the full system, the quantum 
number Q could represent particle number. 

For matrix product states, the introduction of Abelian symmetries has the 
consequence that the A-matrix may be written as (^q^q^)!'^/?,.- Here Qo-, Qh Qr 
are the quantum numbers attached to the local, left effective and right effective basis, 
respectively. The index ai distinguishes different states \Qi,ai) characterized by the 
same quantum number Qi, and similarly for \Qr,Pr) and \Qa,1a)- If A describes, for 
example, the mapping of the |/)-basis of the left block together with the local basis to 
a combined (truncated) |r)-basis, then the only non-zero blocks of the A-matrix are 
those for which Qa + Qi = Qr- For the current site, the total symmetry Qtot of the full 
quantum many-body state manifests itself in that the corresponding A-matrix fulfills 

Ql + Qr + Qa = Qtot- 

For the handling of matrix product states quantum numbers imply a significant 
amount of bookkeeping, i.e. for every coefficient block we have to store its quantum 
number. The benefit is that we can deal with large effective state spaces at reasonable 
numerical cost. The Lanczos algorithm, in particular, takes advantage of the block 
structure. 

Of course, the treatment of Abelian symmetries is generic and not limited to 
only one symmetry. We may incorporate as many symmetries as exist for a given 
Hamiltonian, by writing Q as a vector of the corresponding quantum numbers. 



A. 5. Additional details 



A. 5.1. Derivation of the orthonormality condition The orthonormality condition (A. 7) 



is easily derived by induction. The starting point is condition (A. 6) and we limit to the 
derivation for the left basis. The derivation for the right basis is analogous. 

The induction argument can be initialized with site k = 1 because its effective left 
basis is already orthonormal as it consists only of the vacuum state. Now, consider the 
case that site k has an orthonormal effective left basis and construct the condition for 
site + 1 to have an orthonormal effective left basis: 



fc+1 
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^^Mt^Mj _ (A.38) 



Condition (A. 7) follows with = 6i'^^^i^^^. 



A. 5. 2. Singular value decomposition The singular value decomposition can be seen as 
a generalization of the spectral theorem, i.e. of the eigenvalue decomposition. It is valid 
for any real or complex m x n rectangular matrix. Let M be such a matrix, then it can 
be written in a singular value decomposition 

M = USV^ , (A.39) 

where U is amxm unitary matrix, S am xn matrix with real, nonnegative entries on 
the diagonal and zeros off the diagonal, and V a n x n unitary matrix. The numbers on 
the diagonal of S are called singular values, and there are p = min (n, m) of them. The 
singular values are unique, but U and V are not, in general. It is convenient to truncate 
and reorder these matrices in such a fashion that their dimension are mx p ioi U, px p 
for S (with the singular values ordered in a non-increasing fashion) and nx p for V (i.e. 
p X n for V''). A consequence of this truncation is that f/ or is no longer quadratic and 
unitarity is not defined for such matrices. This property is replaced by column unitarity 
(orthonormal columns) of U and row unitarity (orthonormal rows) for V"^ - no matter 
which one is no longer quadratic. In this article all singular value decompositions are 
understood to be ordered in this fashion. 



A. 5. 3. Numerical costs of index contractions The numerical costs of matrix 
multiplications and index contractions of multi-index objects depend on the dimension 
of both the resulting object and of the contracted indices. In the case of matrix 
multiplications this is quite simple. Consider a n x m matrix Mi multiplied by a m x p 
matrix M2. The result is a n x p matrix M: 

m 

M^, =Y.(Mi\kiM2\, ■ (A.40) 

k=l 

Evidently, each of the n* p matrix elements Mij requires a sum over m products of the 
form {Mi)-f^ {M2)i^j. Thus the process for calculating M1M2 is of order O (nmp). 

The numerical costs of multi-index objects are obtained analogously. Consider two 
multi-index objects. Mi with indices ii, . . . ,in and dimensions piX . . . xpn and M2 with 
indices ji, . . . ,jm and dimensions gi x . . . x q^. If we contract the indices ii and ^2 of 
Ml with the indices ji and j2 of M2 (assuming that pi = qi and p2 = ^2), we obtain the 
multi-index object M: 

Pl P2 
k=l 1=1 

Thus for every entry of M, pi times p2 multiplications have to be done, so that the 
process is of order O {{ps . . . p„) {P1P2) (qs ■ ■ ■ Qm))- 
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